READ ME

This R code is used to take raw respirometry slopes measured on Galaxias maculatus at Deakin’s Marine Research Centre in Queenscliff, Victoria, Australia and calculate their standard, routine, and maximum metabolic rate. These metabolic measurements are made repeatedly on known individuals held under one of two temperature treatments (18°C or 23°C), and paired with repeated mass and length measurements over five month of data collection (May - September, 2023). After reading in all data, a multivariate mixed effects modeling approach is undertaken in order to assess covariance between growth and metabolism traits among- and within-individuals.

LOAD REQUIRED PACKAGES

# installs and loads packages, need to install "pacman" first

pacman::p_load("brms",
               "broom",
               "broom.mixed",
               "car",
               "coda",
               "dplyr",
               "ggeffects",
               "ggiraphExtra",
               "ggplot2",
               "ggpubr",
               "hms",
               "lme4",
               "lmerTest",
               "lubridate",
               "nlme",
               "posterior",
               "readxl",
               "tidyverse"
)

Working directories

wd <- getwd() # getwd tells us what the current wd is, we are using this to drop it in a variable called wd

raw_data_wd <- paste0(wd, "./00_raw_data")  # creates a variable with the name of the wd where raw data are stored

    MRcalcs_wd <- paste0(raw_data_wd, "./MRcalcs")
    
    MR_slopes_wd <- paste0(raw_data_wd, "./MR_Slopes")
    
    MMR_wd <- paste0(raw_data_wd, "./MMR")

metadata_wd <- paste0(wd, "./01_metadata")  # creates a variable with the name of the wd where metadata are stored

figures_wd <- paste0(wd, "./03_figures")  # creates a variable with the name of the wd where figures are to be stored

Read in data, calculate metabolic rates, merge into master data table

#calcSMR function adapted from Chabot, 2016 (https://doi.org/10.1111/jfb.12845)
#note that we have not used the mlnd function as the mean of the lowest 10th percentile was better suited to our data

calcSMR <- function(Y) {
  u <- sort(Y)
  tenpc <- round(0.1 * length(u))
  SD10pc <- sd(u[1:tenpc])
  low10pc = mean(u[(which((u > (mean(u[1:tenpc])-SD10pc)))):(tenpc+which((u > (mean(u[1:tenpc])-SD10pc -u[1]))))])
  return(list(low10pc = low10pc))
}

# Define a function to extract date strings from file names
extract_date <- function(files) {
  gsub(".*?([0-9]{8}).*", "\\1", basename(files))
}

# Specify the directories where your files are located
csv_folder <- MR_slopes_wd
xlsx_folder <- MMR_wd

# Get the list of files in each folder
csv_files <- list.files(csv_folder, pattern = "_MR_slopes.csv", full.names = TRUE)
xlsx_files <- list.files(xlsx_folder, pattern = "_MMR.xlsx", full.names = TRUE)


# Extract date strings from file names
csv_dates <- extract_date(csv_files)
xlsx_dates <- extract_date(xlsx_files)


matching_files <- Map(function(date) {
  csv_file <- csv_files[grep(date, csv_dates)]
  xlsx_file <- xlsx_files[grep(date, xlsx_dates)]
  return(list(csv = csv_file, xlsx = xlsx_file))
}, unique(c(csv_dates, xlsx_dates)))


# Loop through matching files and perform your data processing
for (i in seq_along(matching_files)) {
  date <- unique(csv_dates)[i]
  current_files <- matching_files[[i]]  # Renamed 'files' to 'current_files'
  
  # Process CSV data
  tb_respirometry <- read_csv(current_files$csv) %>%
    rename(
      chamber_ch = Chamber.No,
      ID_fish    = Ind,
      mass       = Mass,
      length     = Length,
      volume_ch  = Ch.Volume,
      DOunit     = DO.unit,
      dateTime   = Date.Time,
      Temp_class = Temp.class,
      slope_wBR  = Slope.with.BR,
      BRSlope    = BR.Slope
    )  %>%
    mutate(
      dateTime       = as.POSIXct(ymd_hms(dateTime)),
      Time           = as_hms(ymd_hms(dateTime)),
      Date           = as.Date(dateTime, format = "%Y/%m/%d"),
    )
  
  
  tb_respirometry <- tb_respirometry %>%
    mutate(
      volume_net = volume_ch - mass,
      MR_wBR     = abs(slope_wBR)*(volume_net/1000)*60*60, #uncomment for mgO2/hr instead
      BR         = BRSlope*(volume_ch/1000)*60*60,
      MR         = MR_wBR + BR,
      BR_perc    = (BR/MR_wBR)*100
    ) 
  
  #Calculate RMR and variance in MR for each fish, create a new table for individual RMR calcs 
  
  tb_rmr <- tb_respirometry %>%
    group_by(chamber_ch) %>%
    arrange(chamber_ch, dateTime) %>%
    slice(3:(n() - 1)) %>%
    ungroup()  %>%
    group_by(
      ID_fish, mass, length, volume_net, chamber_ch, Temp_class
    ) %>% 
    arrange(ID_fish) %>%
    drop_na() %>%
    summarise(
      RMR        = mean(MR),
      RMR_var    = var(MR),
      RMR_perc   = (RMR_var/RMR)*100,
      time_start = dateTime %>% min(),
      time_end   = dateTime %>% max()
    ) 
  
  #Calculate SMR for each fish, create a new table for individual SMR calcs 
  
  tb_smr <-
    tb_respirometry %>%
    group_by(
      ID_fish, mass,length, volume_net, chamber_ch
    ) %>% 
    arrange(ID_fish) %>%
    drop_na() %>%
    summarise(
      SMR        = calcSMR(MR)$low10pc %>% unname(),
      time_start = dateTime %>% min(),
      time_end   = dateTime %>% max()
    ) 
  
    # Process XLSX data
    tb_mmr <- read_excel(current_files$xlsx) %>%
      rename(
        chamber_ch = Chamber.No,
        ID_fish    = Ind,
        mass       = Mass,
        volume_ch  = Ch.Volume,
        DOunit     = DO.unit,
        dateTime   = Date.Time,
        Temp_class = Temp.class,
        slope_wBR  = Slope.with.BR,
        BRSlope    = BR.Slope
      )  %>%
      mutate(
        dateTime       = as.POSIXct(ymd_hms(dateTime)),
        Time           = as_hms(ymd_hms(dateTime)),
        Date           = as.Date(dateTime, format = "%Y/%m/%d"),
      ) %>%
      drop_na()
    
    tb_mmr <- tb_mmr %>%
      mutate(
        volume_net = volume_ch - mass,
        MMR_wBR    = abs(slope_wBR)*(volume_net/1000)*60*60, #uncomment for mgO2/hr instead
        BR         = abs(BRSlope)*(volume_ch/1000)*60*60,
        MMR        = MMR_wBR - BR,
        BR_perc    = (BR/MMR_wBR)*100
      ) %>%
      arrange(ID_fish)
      
    # Combine the processed data
    # Combine the processed data, selecting only specific columns
    tb_MR_master <- tb_rmr %>%
      left_join(select(tb_smr, ID_fish, SMR), by = "ID_fish") %>%
      left_join(select(tb_mmr, ID_fish, MMR), by = "ID_fish") %>%
      select(-matches(".x$")) %>%
      rename_with(~gsub("\\.y$", "", .), matches(".y$"))
    
  # Export data
  setwd(MRcalcs_wd)
  write.csv(tb_MR_master, file = paste0(date, "_MRcalcs.csv"), col.names = NA, row.names = FALSE)
}

# read all "_MRcalcs.csv" files, bind into one table, and add columns with logged values 
file.list <- list.files(MRcalcs_wd)

setwd(MRcalcs_wd)
calcs_all <- lapply(file.list, read_csv)

tb_MRcalcs <- bind_rows(calcs_all, .id="Resp_Day") %>%
  mutate(
    log_SMR_low10pc = log10(SMR),
    log_mass = log10(mass),
    log_RMR = log10(RMR),
    log_MMR = log10(MMR),
    Resp_Day = as.numeric(Resp_Day)
  )

# assign Month category to all timepoints 
tb_MRcalcs$Month <- cut(
  tb_MRcalcs$Resp_Day,
  breaks = c(0, 13, seq(23, max(tb_MRcalcs$Resp_Day) + 10, by = 10)),
  right = FALSE,
  labels = c("April", "May", "June", "July", "August", "September")
)

setwd(raw_data_wd)
biometrics <- read.csv("FinalBiometrics.csv")

ds <- left_join(tb_MRcalcs, biometrics, by = c("ID_fish")) #a table with all resp data and metadata collected in dissection

ds1 <- ds[which(ds$Month != "April"),] # remove April data from analysis due to temperature inconsistency

Multivariate mixed-effects models to test the relationships between growth and metabolism in Galaxias maculatus.

ds1 <- ds1 %>%
  # Filter to keep only ID_fish where a mass exists for Month == "May"
  filter(ID_fish %in% ID_fish[Month == "May" & !is.na(mass)])


month_order <- c("May", "June", "July", "August", "September") # assign chronological order to Month (otherwise alphabetical)
ds1$Month <- factor(ds1$Month, levels = month_order) # assign month_order to factor variable "Month" in ds1

month_order_pop <- c("April", "May", "June", "July", "August", "September") # assign chronological order to Month (otherwise alphabetical)

min_mass = min(ds1$mass)
min_log_mass = min(ds1$log_mass)

min_finalmass = min(((ds1[which(!is.na(ds1$Final_Mass_g)),]$Final_Mass_g)))
min_log_finalmass = min((log10(ds1[which(!is.na(ds1$Final_Mass_g)),]$Final_Mass_g)))


# format data columns and generate centred and log-scaled variables
ds1 <- ds1 %>%
  mutate(
    mass_centre = mass - min_mass,                                  # centred mass  
    log_mass_centre = log_mass - log10(min_mass),                   # log-scaled centred mass  
    Temp_class = as.factor(Temp_class),                             # set temperature treatments as factors
    SMR1 = log10(SMR),                                              # log-scaled SMR
    RMR1 = log10(RMR),                                              # log-scaled RMR
    MMR1 = log10(MMR),                                              # log-scaled MMR
    mass1 = log_mass,                                               # log-scaled mass
    Month_int = case_when(Month == "May" ~ 1,
                          Month == "June" ~ 2,
                          Month == "July" ~ 3,
                          Month == "August" ~ 4,
                          Month == "September" ~ 5),
    Month_centred = Month_int - 1,                                  # left-centred mass so May is the intercept
    Month_centred_factor = as.factor(Month_centred),                # set centred month as factors
    Rack = as.integer(Rack),                                        # set rack ID to integer (for now)
    Bin = as.integer(Bin),                                          # set bin ID to integer (for now)
    finalmass_centre = Final_Mass_g - min_finalmass,                # centre last mass measurement taken
    log_finalmass_centre = log10(Final_Mass_g) - min_log_finalmass, # log-scaled and centred final mass
    log_fatmass = log10(Fat_mass_g)                                 # log-scaled fat mass
  )

ds1 <- tidyr::unite(data = ds1, Rack, Bin, col = "tankID", sep = "_") # merge rack and bin data into one identifier, call it tankID

#read in density data table, call it BinCounts
setwd(raw_data_wd)
BinCounts <- read_excel("BinCounts.xlsx", sheet = "data")

tb_density <- BinCounts %>% 
  group_by_(
    "Month",
    "Rack") %>% 
  summarise(
    Population = sum(Population_tagged)) %>% 
  mutate(
    Month = factor(Month, levels = month_order_pop))

ggplot(BinCounts, aes(x = factor(Month, levels = month_order_pop), y = Population_tagged, group = tankID, color = factor(tankID))) +
  geom_line() +
  geom_point() +
  geom_jitter(height = 0.1, width = 0.25) +
  facet_grid(.~Rack) +
  theme_classic()

#merge BinCounts and ds1
ds1 <- left_join(ds1, BinCounts %>% dplyr::select("Month", "tankID", "Population"), by = c("Month", "tankID"))

#reformat Sex variable so that both unknown and immature fish are "NA" 
ds1$Sex <- ifelse(ds1$Sex == "I", NA, ds1$Sex)

# rename "population" as density, set tankID as factor
ds1 <- ds1 %>%
  rename(
    Density = Population
  ) %>%
  mutate(
    tankID = as.factor(tankID)
  )

# check structure of ds1
str(ds1)
## tibble [624 × 41] (S3: tbl_df/tbl/data.frame)
##  $ Resp_Day            : num [1:624] 13 13 13 13 13 13 13 13 13 13 ...
##  $ ID_fish             : num [1:624] 81 82 83 84 85 86 87 88 89 90 ...
##  $ chamber_ch          : chr [1:624] "A1" "B3" "B4" "A2" ...
##  $ Temp_class          : Factor w/ 2 levels "18","23": 2 2 2 2 2 2 2 2 2 2 ...
##  $ RMR                 : num [1:624] 0.952 0.412 0.712 1.016 0.876 ...
##  $ RMR_var             : num [1:624] 0.0149 0.0135 0.0429 0.0696 0.0454 ...
##  $ RMR_perc            : num [1:624] 1.57 3.28 6.02 6.85 5.18 ...
##  $ time_start          : POSIXct[1:624], format: "2023-05-22 17:53:00" "2023-05-22 17:53:00" ...
##  $ time_end            : POSIXct[1:624], format: "2023-05-23 08:03:00" "2023-05-23 08:03:00" ...
##  $ mass                : num [1:624] 2.65 0.81 1.78 2.73 2.29 0.98 2.04 2.43 2.31 1.65 ...
##  $ length              : num [1:624] 77 62 73 78 80 64 80 77 78 68 ...
##  $ volume_net          : num [1:624] 707 709 708 707 708 ...
##  $ SMR                 : num [1:624] 0.845 0.284 0.53 0.804 0.707 ...
##  $ MMR                 : num [1:624] 1.43 1.11 1.44 1.64 1.69 ...
##  $ log_SMR_low10pc     : num [1:624] -0.0733 -0.5461 -0.2758 -0.0948 -0.1506 ...
##  $ log_mass            : num [1:624] 0.4232 -0.0915 0.2504 0.4362 0.3598 ...
##  $ log_RMR             : num [1:624] -0.02124 -0.38552 -0.14722 0.00674 -0.05726 ...
##  $ log_MMR             : num [1:624] 0.1544 0.0458 0.1591 0.2154 0.2275 ...
##  $ Month               : chr [1:624] "May" "May" "May" "May" ...
##  $ tankID              : Factor w/ 20 levels "1_10","1_11",..: 17 17 17 17 17 17 17 17 18 18 ...
##  $ TagID               : chr [1:624] "YL" "YR" "BrL" "BrR" ...
##  $ Final_Mass_g        : num [1:624] 4.62 3.65 2.64 6.33 5.19 ...
##  $ Final_Length_mm     : int [1:624] 94 86 84 98 100 95 90 89 91 98 ...
##  $ Gonad_mass_g        : num [1:624] 0.035 0.022 0.018 0.039 0.052 0.028 0.362 0.033 0.015 0.093 ...
##  $ Sex                 : chr [1:624] "F" "M" "F" "M" ...
##  $ Fat_mass_g          : num [1:624] 0.315 0.123 0.062 0.314 0.214 0.091 0.076 0.044 0.15 0.209 ...
##  $ Ventricle_mass_g    : chr [1:624] "0.017" "" "" "" ...
##  $ Otolith_age_y       : num [1:624] 1 NA NA NA NA 1 1 NA NA 1 ...
##  $ mass_centre         : num [1:624] 2.07 0.23 1.2 2.15 1.71 0.4 1.46 1.85 1.73 1.07 ...
##  $ log_mass_centre     : num [1:624] 0.66 0.145 0.487 0.673 0.596 ...
##  $ SMR1                : num [1:624] -0.0733 -0.5461 -0.2758 -0.0948 -0.1506 ...
##  $ RMR1                : num [1:624] -0.02124 -0.38552 -0.14722 0.00674 -0.05726 ...
##  $ MMR1                : num [1:624] 0.1544 0.0458 0.1591 0.2154 0.2275 ...
##  $ mass1               : num [1:624] 0.4232 -0.0915 0.2504 0.4362 0.3598 ...
##  $ Month_int           : num [1:624] 1 1 1 1 1 1 1 1 1 1 ...
##  $ Month_centred       : num [1:624] 0 0 0 0 0 0 0 0 0 0 ...
##  $ Month_centred_factor: Factor w/ 5 levels "0","1","2","3",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ finalmass_centre    : num [1:624] 3.3 2.33 1.32 5.01 3.87 ...
##  $ log_finalmass_centre: num [1:624] 0.544 0.442 0.301 0.681 0.594 ...
##  $ log_fatmass         : num [1:624] -0.502 -0.91 -1.208 -0.503 -0.67 ...
##  $ Density             : num [1:624] 8 8 8 8 8 8 8 8 5 5 ...
head(ds1)
## # A tibble: 6 × 41
##   Resp_Day ID_fish chamber_ch Temp_class   RMR RMR_var RMR_perc
##      <dbl>   <dbl> <chr>      <fct>      <dbl>   <dbl>    <dbl>
## 1       13      81 A1         23         0.952  0.0149     1.57
## 2       13      82 B3         23         0.412  0.0135     3.28
## 3       13      83 B4         23         0.712  0.0429     6.02
## 4       13      84 A2         23         1.02   0.0696     6.85
## 5       13      85 A4         23         0.876  0.0454     5.18
## 6       13      86 A3         23         0.751  0.0504     6.71
## # ℹ 34 more variables: time_start <dttm>, time_end <dttm>, mass <dbl>,
## #   length <dbl>, volume_net <dbl>, SMR <dbl>, MMR <dbl>,
## #   log_SMR_low10pc <dbl>, log_mass <dbl>, log_RMR <dbl>, log_MMR <dbl>,
## #   Month <chr>, tankID <fct>, TagID <chr>, Final_Mass_g <dbl>,
## #   Final_Length_mm <int>, Gonad_mass_g <dbl>, Sex <chr>, Fat_mass_g <dbl>,
## #   Ventricle_mass_g <chr>, Otolith_age_y <dbl>, mass_centre <dbl>,
## #   log_mass_centre <dbl>, SMR1 <dbl>, RMR1 <dbl>, MMR1 <dbl>, mass1 <dbl>, …
# check sex ratio of population
tb_sexratio <- ds1[which(ds1$Month == "September"),] %>%
  group_by_(
    "Temp_class",
    "Sex"
  ) %>% 
  count() 

sexorder <- c("Immature", "Female", "Male")

tb_sexratio$Sex <- as.character(tb_sexratio$Sex)
tb_sexratio$Sex[tb_sexratio$Sex==""] <- NA
tb_sexratio$Sex[is.na(tb_sexratio$Sex)] <- "Immature"
tb_sexratio$Sex[tb_sexratio$Sex == "F"] <- "Female"
tb_sexratio$Sex[tb_sexratio$Sex == "M"] <- "Male"
tb_sexratio <- tb_sexratio %>%
  mutate(
    Sex = factor(Sex, levels = sexorder))

fig3b1 <- ggplot(tb_sexratio, aes(x = " ", y = n,  fill = Sex)) +
  geom_bar(position = "stack", stat = "identity") +
  scale_fill_manual(values = c("gold3", "darkgreen", "purple")) +
  labs(x = "Whole Population", y = "Count (absolute)", fill = "Sex") +
  theme_classic()

fig3b2 <- ggplot(tb_sexratio, aes(x = factor(Temp_class), y = n, fill = Sex)) +
  geom_bar(position = "fill", stat = "identity") +
  scale_fill_manual(values = c("gold3", "darkgreen", "purple")) +
  labs(x = expression("Temperature ("*degree*C*")"), y = "Proportion", fill = "Sex") +
  theme_classic()

fig3b <- ggarrange(print(fig3b1), print(fig3b2), ncol = 2, nrow = 1, widths = c(1, 2), align = "v", common.legend = TRUE, legend = "right")

fig3b

GROWTH ANALYSES

Visualize growth data

hist(ds1$mass) # plot all mass values to see mass range and most frequent values

ds1$Month <- factor(ds1$Month, levels = month_order)              # reassign month_order to factor variable "Month" in ds1

ds1$Sex[ds1$Sex == ""] <- NA
ds1$Sex <- factor(ds1$Sex)


# Plot mass over time, connect by ID_fish to show individual trends

ggplot(ds1, aes(x = Month, y = mass, group = ID_fish)) +          # plot raw data by ID
  geom_point(size = 2) +                                          # size of data points
  geom_line() +                                                   # connect points by individual ID
  facet_wrap(. ~ Temp_class) +                                    # produces one plot per temperature treatment
  theme_bw()

Log-linearize mass data

We can also visualise this by plotting log10 mass over log10_x (note that we cannot use log10_y as some masses were < 1 initially)

ggplot(ds1, aes(x = Month_int, y = log_mass, group = ID_fish)) +  # plot logged mass data by ID
  geom_point(size = 2) +                                          # size of data points
  geom_line() +                                                   # connect points by individual ID
  facet_wrap(. ~ Temp_class) +                                    # produces one plot per temperature treatment
  theme_bw() +
  scale_x_log10()                                                 # display x-axis on log scale

Univariate growth models for two temperature treatments

ds18 <- ds1[which(ds1$Temp_class == "18"),]                   # subset ds for only fish at 18C
ds23 <- ds1[which(ds1$Temp_class == "23"),]                   # subset ds for only fish at 23C

mass18_mod<-lme4::lmer(log_mass ~ Month + (1 + Month_int|ID_fish), # model structure, 18C only
           data=ds18, REML = T)

summary(mass18_mod)                                 # generate model output
## Linear mixed model fit by REML ['lmerMod']
## Formula: log_mass ~ Month + (1 + Month_int | ID_fish)
##    Data: ds18
## 
## REML criterion at convergence: -492
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.8214 -0.4648  0.0139  0.4865  2.3450 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  ID_fish  (Intercept) 0.040770 0.20191       
##           Month_int   0.002038 0.04515  -0.57
##  Residual             0.003528 0.05939       
## Number of obs: 317, groups:  ID_fish, 72
## 
## Fixed effects:
##                Estimate Std. Error t value
## (Intercept)     0.32789    0.02232  14.689
## MonthJune       0.15461    0.01147  13.478
## MonthJuly       0.22113    0.01498  14.761
## MonthAugust     0.15515    0.01989   7.799
## MonthSeptember  0.36370    0.02487  14.622
## 
## Correlation of Fixed Effects:
##             (Intr) MnthJn MnthJl MnthAg
## MonthJune   -0.364                     
## MonthJuly   -0.412  0.629              
## MonthAugust -0.410  0.600  0.770       
## MonthSptmbr -0.408  0.581  0.776  0.865
plot(mass18_mod)                                    # residuals plot

mass23_mod<-lme4::lmer(log_mass ~ Month + (1 + Month_int|ID_fish), # model structure, 23C only
                 data=ds23, REML=T)

summary(mass23_mod)                                 # generate model output
## Linear mixed model fit by REML ['lmerMod']
## Formula: log_mass ~ Month + (1 + Month_int | ID_fish)
##    Data: ds23
## 
## REML criterion at convergence: -515.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.9289 -0.3941 -0.0018  0.5057  4.8187 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  ID_fish  (Intercept) 0.040995 0.20247       
##           Month_int   0.001445 0.03802  -0.46
##  Residual             0.003138 0.05602       
## Number of obs: 307, groups:  ID_fish, 64
## 
## Fixed effects:
##                Estimate Std. Error t value
## (Intercept)     0.21996    0.02451   8.974
## MonthJune       0.12091    0.01123  10.771
## MonthJuly       0.21448    0.01402  15.297
## MonthAugust     0.22778    0.01773  12.848
## MonthSeptember  0.35157    0.02193  16.028
## 
## Correlation of Fixed Effects:
##             (Intr) MnthJn MnthJl MnthAg
## MonthJune   -0.298                     
## MonthJuly   -0.335  0.623              
## MonthAugust -0.341  0.611  0.773       
## MonthSptmbr -0.337  0.589  0.778  0.858
plot(mass23_mod)                                    # residuals plot

Combined growth model

# mass model including both temperature treatments, and their interaction with Month
mass_mod <- lmer(log_mass ~ Month*as.factor(Temp_class) + Sex  + Density + (1 + Month_int|ID_fish), 
                 data=ds1, REML=T)

summary(mass_mod)   # generate model output
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log_mass ~ Month * as.factor(Temp_class) + Sex + Density + (1 +  
##     Month_int | ID_fish)
##    Data: ds1
## 
## REML criterion at convergence: -732.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6965 -0.4297  0.0123  0.4659  4.5929 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  ID_fish  (Intercept) 0.032801 0.18111       
##           Month_int   0.001623 0.04029  -0.54
##  Residual             0.003489 0.05906       
## Number of obs: 453, groups:  ID_fish, 91
## 
## Fixed effects:
##                                          Estimate Std. Error         df t value
## (Intercept)                              0.287001   0.051703 379.511550   5.551
## MonthJune                                0.157704   0.012697 354.616467  12.420
## MonthJuly                                0.228444   0.016168 237.622029  14.130
## MonthAugust                              0.167240   0.021107 155.507265   7.923
## MonthSeptember                           0.379389   0.025536 111.141150  14.857
## as.factor(Temp_class)23                 -0.059851   0.037534  99.223469  -1.595
## SexM                                    -0.073375   0.033162  88.017749  -2.213
## Density                                  0.011296   0.005978 429.998459   1.890
## MonthJune:as.factor(Temp_class)23       -0.024598   0.019947 352.647746  -1.233
## MonthJuly:as.factor(Temp_class)23       -0.010083   0.025097 226.080174  -0.402
## MonthAugust:as.factor(Temp_class)23      0.060749   0.032152 137.635190   1.889
## MonthSeptember:as.factor(Temp_class)23  -0.023724   0.039367 100.653936  -0.603
##                                        Pr(>|t|)    
## (Intercept)                            5.35e-08 ***
## MonthJune                               < 2e-16 ***
## MonthJuly                               < 2e-16 ***
## MonthAugust                            4.15e-13 ***
## MonthSeptember                          < 2e-16 ***
## as.factor(Temp_class)23                  0.1140    
## SexM                                     0.0295 *  
## Density                                  0.0595 .  
## MonthJune:as.factor(Temp_class)23        0.2183    
## MonthJuly:as.factor(Temp_class)23        0.6882    
## MonthAugust:as.factor(Temp_class)23      0.0609 .  
## MonthSeptember:as.factor(Temp_class)23   0.5481    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##                (Intr) MnthJn MnthJl MnthAg MnthSp a.(T_) SexM   Densty
## MonthJune      -0.281                                                 
## MonthJuly      -0.384  0.632                                          
## MonthAugust    -0.473  0.615  0.787                                   
## MonthSptmbr    -0.423  0.599  0.794  0.873                            
## as.fc(T_)23    -0.358  0.239  0.273  0.282  0.278                     
## SexM           -0.215 -0.004 -0.006 -0.007 -0.006 -0.098              
## Density        -0.861  0.139  0.239  0.345  0.285  0.116 -0.021       
## MnthJn:.(T_)23  0.143 -0.631 -0.392 -0.377 -0.370 -0.364  0.002 -0.047
## MnthJl:.(T_)23  0.219 -0.402 -0.636 -0.496 -0.502 -0.414  0.003 -0.121
## MnA:.(T_)23     0.288 -0.400 -0.511 -0.648 -0.566 -0.425  0.004 -0.201
## MnS:.(T_)23     0.245 -0.384 -0.507 -0.555 -0.639 -0.419  0.003 -0.151
##                MnthJn:.(T_)23 MnthJl:.(T_)23 MA:.(T
## MonthJune                                          
## MonthJuly                                          
## MonthAugust                                        
## MonthSptmbr                                        
## as.fc(T_)23                                        
## SexM                                               
## Density                                            
## MnthJn:.(T_)23                                     
## MnthJl:.(T_)23  0.625                              
## MnA:.(T_)23     0.608          0.778               
## MnS:.(T_)23     0.591          0.785          0.864
plot(mass_mod)      # residuals plot

ranova(mass_mod)    # ANOVA table output to test random effects
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log_mass ~ Month + as.factor(Temp_class) + Sex + Density + (1 + Month_int | ID_fish) + Month:as.factor(Temp_class)
##                                        npar logLik     AIC    LRT Df Pr(>Chisq)
## <none>                                   16 366.13 -700.25                     
## Month_int in (1 + Month_int | ID_fish)   14 305.38 -582.77 121.48  2  < 2.2e-16
##                                           
## <none>                                    
## Month_int in (1 + Month_int | ID_fish) ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# mass model including both temperature treatments, and their interaction with Month
mass_mod_noMonth <- lmer(log_mass ~ Month*as.factor(Temp_class) + Sex + Density + (1|ID_fish), 
                 data=ds1, REML=T)

summary(mass_mod_noMonth)   # generate model output
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log_mass ~ Month * as.factor(Temp_class) + Sex + Density + (1 |  
##     ID_fish)
##    Data: ds1
## 
## REML criterion at convergence: -610.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.6634 -0.5041  0.0313  0.4740  3.9635 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID_fish  (Intercept) 0.022746 0.15082 
##  Residual             0.007567 0.08699 
## Number of obs: 453, groups:  ID_fish, 91
## 
## Fixed effects:
##                                          Estimate Std. Error         df t value
## (Intercept)                              0.304207   0.053767 382.535854   5.658
## MonthJune                                0.157025   0.016791 355.224021   9.352
## MonthJuly                                0.227033   0.017176 361.268479  13.218
## MonthAugust                              0.164380   0.018271 377.508492   8.997
## MonthSeptember                           0.376529   0.018271 377.508492  20.608
## as.factor(Temp_class)23                 -0.061639   0.037779 137.305730  -1.632
## SexM                                    -0.072332   0.033449  88.014836  -2.162
## Density                                  0.008948   0.006287 440.998220   1.423
## MonthJune:as.factor(Temp_class)23       -0.024244   0.026455 353.637257  -0.916
## MonthJuly:as.factor(Temp_class)23       -0.008998   0.026629 355.432235  -0.338
## MonthAugust:as.factor(Temp_class)23      0.063283   0.027233 362.066784   2.324
## MonthSeptember:as.factor(Temp_class)23  -0.021386   0.027107 360.854223  -0.789
##                                        Pr(>|t|)    
## (Intercept)                            3.01e-08 ***
## MonthJune                               < 2e-16 ***
## MonthJuly                               < 2e-16 ***
## MonthAugust                             < 2e-16 ***
## MonthSeptember                          < 2e-16 ***
## as.factor(Temp_class)23                  0.1051    
## SexM                                     0.0333 *  
## Density                                  0.1554    
## MonthJune:as.factor(Temp_class)23        0.3601    
## MonthJuly:as.factor(Temp_class)23        0.7356    
## MonthAugust:as.factor(Temp_class)23      0.0207 *  
## MonthSeptember:as.factor(Temp_class)23   0.4307    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##                (Intr) MnthJn MnthJl MnthAg MnthSp a.(T_) SexM   Densty
## MonthJune      -0.248                                                 
## MonthJuly      -0.355  0.506                                          
## MonthAugust    -0.505  0.495  0.538                                   
## MonthSptmbr    -0.505  0.495  0.538  0.588                            
## as.fc(T_)23    -0.355  0.230  0.241  0.250  0.250                     
## SexM           -0.206 -0.004 -0.007 -0.010 -0.010 -0.099              
## Density        -0.871  0.110  0.237  0.419  0.419  0.121 -0.023       
## MnthJn:.(T_)23  0.129 -0.631 -0.313 -0.300 -0.300 -0.353  0.002 -0.037
## MnthJl:.(T_)23  0.200 -0.323 -0.637 -0.333 -0.333 -0.360  0.004 -0.120
## MnA:.(T_)23     0.311 -0.328 -0.353 -0.657 -0.381 -0.368  0.006 -0.249
## MnS:.(T_)23     0.296 -0.328 -0.350 -0.375 -0.652 -0.368  0.005 -0.231
##                MnthJn:.(T_)23 MnthJl:.(T_)23 MA:.(T
## MonthJune                                          
## MonthJuly                                          
## MonthAugust                                        
## MonthSptmbr                                        
## as.fc(T_)23                                        
## SexM                                               
## Density                                            
## MnthJn:.(T_)23                                     
## MnthJl:.(T_)23  0.499                              
## MnA:.(T_)23     0.492          0.509               
## MnS:.(T_)23     0.494          0.510          0.529
plot(mass_mod_noMonth)      # residuals plot

ranova(mass_mod_noMonth)    # ANOVA table output to test random effects
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log_mass ~ Month + as.factor(Temp_class) + Sex + Density + (1 | ID_fish) + Month:as.factor(Temp_class)
##               npar logLik     AIC    LRT Df Pr(>Chisq)    
## <none>          14 305.38 -582.77                         
## (1 | ID_fish)   13 121.89 -217.78 366.99  1  < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# mass model with random intercepts by temp treatment
mass_mod_trtint <- lmer(log_mass ~ Month*as.factor(Temp_class) + Sex + Density + (0 + Month_int|Temp_class), 
                 data=ds1, REML=T)

summary(mass_mod_trtint)   # generate model output
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log_mass ~ Month * as.factor(Temp_class) + Sex + Density + (0 +  
##     Month_int | Temp_class)
##    Data: ds1
## 
## REML criterion at convergence: -243.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.86639 -0.70696 -0.01277  0.67903  2.94871 
## 
## Random effects:
##  Groups     Name      Variance Std.Dev.
##  Temp_class Month_int 0.03024  0.1739  
##  Residual             0.03013  0.1736  
## Number of obs: 453, groups:  Temp_class, 2
## 
## Fixed effects:
##                                          Estimate Std. Error         df t value
## (Intercept)                              0.301995   0.181676 441.000000   1.662
## MonthJune                                0.159454   0.177058 441.000000   0.901
## MonthJuly                                0.229581   0.349405 441.000000   0.657
## MonthAugust                              0.164793   0.522797 441.000000   0.315
## MonthSeptember                           0.376942   0.696419 441.000000   0.541
## as.factor(Temp_class)23                 -0.061299   0.248776 441.000000  -0.246
## SexM                                    -0.073232   0.016703 441.000000  -4.384
## Density                                  0.009287   0.006273 441.000000   1.480
## MonthJune:as.factor(Temp_class)23       -0.026626   0.251518 441.000000  -0.106
## MonthJuly:as.factor(Temp_class)23       -0.011500   0.494684 441.000000  -0.023
## MonthAugust:as.factor(Temp_class)23      0.062917   0.739690 441.000000   0.085
## MonthSeptember:as.factor(Temp_class)23  -0.021724   0.985139 441.000000  -0.022
##                                        Pr(>|t|)    
## (Intercept)                              0.0972 .  
## MonthJune                                0.3683    
## MonthJuly                                0.5115    
## MonthAugust                              0.7527    
## MonthSeptember                           0.5886    
## as.factor(Temp_class)23                  0.8055    
## SexM                                   1.46e-05 ***
## Density                                  0.1395    
## MonthJune:as.factor(Temp_class)23        0.9157    
## MonthJuly:as.factor(Temp_class)23        0.9815    
## MonthAugust:as.factor(Temp_class)23      0.9323    
## MonthSeptember:as.factor(Temp_class)23   0.9824    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##                (Intr) MnthJn MnthJl MnthAg MnthSp a.(T_) SexM   Densty
## MonthJune       0.920                                                 
## MonthJuly       0.941  0.987                                          
## MonthAugust     0.946  0.986  0.996                                   
## MonthSptmbr     0.949  0.986  0.997  0.998                            
## as.fc(T_)23    -0.686 -0.674 -0.689 -0.693 -0.695                     
## SexM           -0.021 -0.001 -0.001 -0.001 -0.001 -0.008              
## Density        -0.257  0.011  0.012  0.015  0.011  0.018 -0.047       
## MnthJn:.(T_)23 -0.649 -0.704 -0.694 -0.694 -0.694  0.944  0.001 -0.004
## MnthJl:.(T_)23 -0.665 -0.697 -0.706 -0.704 -0.704  0.972  0.001 -0.007
## MnA:.(T_)23    -0.669 -0.697 -0.704 -0.707 -0.706  0.978  0.000 -0.009
## MnS:.(T_)23    -0.671 -0.697 -0.704 -0.706 -0.707  0.981  0.000 -0.006
##                MnthJn:.(T_)23 MnthJl:.(T_)23 MA:.(T
## MonthJune                                          
## MonthJuly                                          
## MonthAugust                                        
## MonthSptmbr                                        
## as.fc(T_)23                                        
## SexM                                               
## Density                                            
## MnthJn:.(T_)23                                     
## MnthJl:.(T_)23  0.983                              
## MnA:.(T_)23     0.983          0.996               
## MnS:.(T_)23     0.982          0.996          0.998
plot(mass_mod_trtint)      # residuals plot

ranova(mass_mod_trtint)    # ANOVA table output to test random effects
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log_mass ~ Month + as.factor(Temp_class) + Sex + Density + (0 + Month_int | Temp_class) + Month:as.factor(Temp_class)
##                                           npar logLik     AIC         LRT Df
## <none>                                      14 121.89 -215.78               
## Month_int in (0 + Month_int | Temp_class)   14 121.89 -215.78 -1.7053e-13  0
##                                           Pr(>Chisq)
## <none>                                              
## Month_int in (0 + Month_int | Temp_class)
mass_mod_sex <- lmer(log_mass ~ 0 + as.factor(Sex):as.factor(Temp_class) + as.factor(Sex):as.factor(Temp_class):Month + Density + (1 + Month_int|ID_fish), 
                 data=ds1, REML=T)
summary(mass_mod_sex)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## log_mass ~ 0 + as.factor(Sex):as.factor(Temp_class) + as.factor(Sex):as.factor(Temp_class):Month +  
##     Density + (1 + Month_int | ID_fish)
##    Data: ds1
## 
## REML criterion at convergence: -706.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.4000 -0.4181  0.0150  0.4378  4.4184 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev. Corr 
##  ID_fish  (Intercept) 0.033320 0.18254       
##           Month_int   0.001501 0.03874  -0.56
##  Residual             0.003432 0.05858       
## Number of obs: 453, groups:  ID_fish, 91
## 
## Fixed effects:
##                                                         Estimate Std. Error
## Density                                                1.228e-02  5.915e-03
## as.factor(Sex)F:as.factor(Temp_class)18                2.817e-01  5.316e-02
## as.factor(Sex)M:as.factor(Temp_class)18                2.026e-01  5.901e-02
## as.factor(Sex)F:as.factor(Temp_class)23                2.237e-01  5.597e-02
## as.factor(Sex)M:as.factor(Temp_class)23                1.435e-01  5.888e-02
## as.factor(Sex)F:as.factor(Temp_class)18:MonthJune      1.598e-01  1.571e-02
## as.factor(Sex)M:as.factor(Temp_class)18:MonthJune      1.544e-01  2.050e-02
## as.factor(Sex)F:as.factor(Temp_class)23:MonthJune      1.078e-01  2.098e-02
## as.factor(Sex)M:as.factor(Temp_class)23:MonthJune      1.617e-01  2.225e-02
## as.factor(Sex)F:as.factor(Temp_class)18:MonthJuly      2.577e-01  1.971e-02
## as.factor(Sex)M:as.factor(Temp_class)18:MonthJuly      1.797e-01  2.557e-02
## as.factor(Sex)F:as.factor(Temp_class)23:MonthJuly      1.689e-01  2.602e-02
## as.factor(Sex)M:as.factor(Temp_class)23:MonthJuly      2.739e-01  2.757e-02
## as.factor(Sex)F:as.factor(Temp_class)18:MonthAugust    1.978e-01  2.530e-02
## as.factor(Sex)M:as.factor(Temp_class)18:MonthAugust    1.170e-01  3.257e-02
## as.factor(Sex)F:as.factor(Temp_class)23:MonthAugust    1.757e-01  3.274e-02
## as.factor(Sex)M:as.factor(Temp_class)23:MonthAugust    2.868e-01  3.466e-02
## as.factor(Sex)F:as.factor(Temp_class)18:MonthSeptember 4.041e-01  3.066e-02
## as.factor(Sex)M:as.factor(Temp_class)18:MonthSeptember 3.394e-01  3.983e-02
## as.factor(Sex)F:as.factor(Temp_class)23:MonthSeptember 2.912e-01  4.032e-02
## as.factor(Sex)M:as.factor(Temp_class)23:MonthSeptember 4.282e-01  4.265e-02
##                                                               df t value
## Density                                                4.209e+02   2.076
## as.factor(Sex)F:as.factor(Temp_class)18                3.488e+02   5.299
## as.factor(Sex)M:as.factor(Temp_class)18                2.778e+02   3.434
## as.factor(Sex)F:as.factor(Temp_class)23                2.403e+02   3.997
## as.factor(Sex)M:as.factor(Temp_class)23                2.379e+02   2.437
## as.factor(Sex)F:as.factor(Temp_class)18:MonthJune      3.451e+02  10.173
## as.factor(Sex)M:as.factor(Temp_class)18:MonthJune      3.442e+02   7.530
## as.factor(Sex)F:as.factor(Temp_class)23:MonthJune      3.436e+02   5.136
## as.factor(Sex)M:as.factor(Temp_class)23:MonthJune      3.443e+02   7.268
## as.factor(Sex)F:as.factor(Temp_class)18:MonthJuly      2.361e+02  13.070
## as.factor(Sex)M:as.factor(Temp_class)18:MonthJuly      2.266e+02   7.028
## as.factor(Sex)F:as.factor(Temp_class)23:MonthJuly      2.229e+02   6.491
## as.factor(Sex)M:as.factor(Temp_class)23:MonthJuly      2.238e+02   9.937
## as.factor(Sex)F:as.factor(Temp_class)18:MonthAugust    1.475e+02   7.821
## as.factor(Sex)M:as.factor(Temp_class)18:MonthAugust    1.366e+02   3.593
## as.factor(Sex)F:as.factor(Temp_class)23:MonthAugust    1.287e+02   5.365
## as.factor(Sex)M:as.factor(Temp_class)23:MonthAugust    1.292e+02   8.273
## as.factor(Sex)F:as.factor(Temp_class)18:MonthSeptember 1.054e+02  13.181
## as.factor(Sex)M:as.factor(Temp_class)18:MonthSeptember 9.942e+01   8.523
## as.factor(Sex)F:as.factor(Temp_class)23:MonthSeptember 9.524e+01   7.221
## as.factor(Sex)M:as.factor(Temp_class)23:MonthSeptember 9.541e+01  10.040
##                                                        Pr(>|t|)    
## Density                                                0.038492 *  
## as.factor(Sex)F:as.factor(Temp_class)18                2.07e-07 ***
## as.factor(Sex)M:as.factor(Temp_class)18                0.000686 ***
## as.factor(Sex)F:as.factor(Temp_class)23                8.53e-05 ***
## as.factor(Sex)M:as.factor(Temp_class)23                0.015541 *  
## as.factor(Sex)F:as.factor(Temp_class)18:MonthJune       < 2e-16 ***
## as.factor(Sex)M:as.factor(Temp_class)18:MonthJune      4.50e-13 ***
## as.factor(Sex)F:as.factor(Temp_class)23:MonthJune      4.70e-07 ***
## as.factor(Sex)M:as.factor(Temp_class)23:MonthJune      2.46e-12 ***
## as.factor(Sex)F:as.factor(Temp_class)18:MonthJuly       < 2e-16 ***
## as.factor(Sex)M:as.factor(Temp_class)18:MonthJuly      2.43e-11 ***
## as.factor(Sex)F:as.factor(Temp_class)23:MonthJuly      5.45e-10 ***
## as.factor(Sex)M:as.factor(Temp_class)23:MonthJuly       < 2e-16 ***
## as.factor(Sex)F:as.factor(Temp_class)18:MonthAugust    9.26e-13 ***
## as.factor(Sex)M:as.factor(Temp_class)18:MonthAugust    0.000455 ***
## as.factor(Sex)F:as.factor(Temp_class)23:MonthAugust    3.63e-07 ***
## as.factor(Sex)M:as.factor(Temp_class)23:MonthAugust    1.39e-13 ***
## as.factor(Sex)F:as.factor(Temp_class)18:MonthSeptember  < 2e-16 ***
## as.factor(Sex)M:as.factor(Temp_class)18:MonthSeptember 1.75e-13 ***
## as.factor(Sex)F:as.factor(Temp_class)23:MonthSeptember 1.26e-10 ***
## as.factor(Sex)M:as.factor(Temp_class)23:MonthSeptember  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
plot(mass_mod_sex)

Plot predicted mean level trend in mass over time, by temperature

pred <- ggpredict(mass_mod, terms = c("Month", "Temp_class")) # generate trend lines for the two temperature treatments

ggplot() +
    geom_errorbar(
      data = pred, 
      aes(x=x, ymin = conf.low, ymax = conf.high, color = group, group = group), width = 0.2) +   # add 95% confidence interval bars
  geom_point(
    data = pred, 
    aes(x = x, y = predicted, color = group, group = group)) +                                    # predicted level by temp*month
  geom_line(
    data = pred, 
    aes(x = x, y = predicted, color = group, group = group)) +                                    # show level trend over time
  labs(
    x = "Month", 
    y = "Predicted log_mass", 
    color = "Temp_class") +
  theme_bw() +
  scale_color_manual(values = c("#78B7C5", "#F21A00"))

F test on interaction effect

anova(mass_mod, type = "3") # type 3 anova because we expect a meaningful interaction
## Type III Analysis of Variance Table with Satterthwaite's method
##                              Sum Sq Mean Sq NumDF  DenDF  F value    Pr(>F)    
## Month                       1.68541 0.42135     4 205.57 120.7818 < 2.2e-16 ***
## as.factor(Temp_class)       0.01090 0.01090     1  88.16   3.1235   0.08063 .  
## Sex                         0.01708 0.01708     1  88.02   4.8958   0.02951 *  
## Density                     0.01246 0.01246     1 430.00   3.5706   0.05948 .  
## Month:as.factor(Temp_class) 0.10282 0.02571     4 201.73   7.3687 1.472e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Generate fitted data

z_mass<-augment(mass_mod)                             # get variables and fitted values from model
head(z_mass)
## # A tibble: 6 × 19
##   .rownames log_mass Month `as.factor(Temp_class)` Sex   Density Month_int
##   <chr>        <dbl> <fct> <fct>                   <fct>   <dbl>     <dbl>
## 1 1           0.423  May   23                      F           8         1
## 2 2          -0.0915 May   23                      M           8         1
## 3 3           0.250  May   23                      F           8         1
## 4 4           0.436  May   23                      M           8         1
## 5 5           0.360  May   23                      M           8         1
## 6 7           0.310  May   23                      M           8         1
## # ℹ 12 more variables: ID_fish <dbl>, .fitted <dbl>, .resid <dbl>, .hat <dbl>,
## #   .cooksd <dbl>, .fixed <dbl>, .mu <dbl>, .offset <dbl>, .sqrtXwt <dbl>,
## #   .sqrtrwt <dbl>, .weights <dbl>, .wtres <dbl>
##### check assumptions: plot all residual values against the fitted 
ggplot(z_mass, aes(x = .fitted, y = .resid)) +   # same plot as before
  geom_point(size = 2) 

FIGURE 1

fig1.1 <- ggplot(z_mass, aes(x = Month, y = log_mass, group = ID_fish)) + 
  geom_line(aes(y = .fitted, color = `as.factor(Temp_class)`),             # note this will look odd if there are additional fixed effects      
            alpha = 0.6) +
  theme_classic() +
  labs(x = "Month", y = expression(log[10]~(mass~(g))), color = expression("Temperature ("*degree*C*")")) +
  ylim(-0.24,1) +
  scale_color_manual(values  = c("#78B7C5", "#F21A00"))

fig1.2 <- ggplot() +
 geom_errorbar(
      data = pred, 
      aes(x=x, ymin = conf.low, ymax = conf.high, color = group, group = group), width = 0.2, alpha = 0.3) +      # add 95% confidence interval bars
  geom_line(
    data = pred, 
    aes(x = x, y = predicted, color = group, group = group), 
            size = 1) +
  labs(
    x = "Month", 
    y = expression(log[10]~(mass~(g))), 
    color = expression("Temperature ("*degree*C*")")) +
  theme_classic() +
  ylim(-0.24,1) +
  scale_color_manual(values = c("#78B7C5", "#F21A00"))

fig1 <- ggarrange(print(fig1.1 + rremove("xlab") + rremove ("ylab")), 
                  print(fig1.2 + rremove("xlab") + rremove ("ylab")), 
                           ncol = 2, nrow = 1, widths = c(1, 1.1),
                           align = "v", 
                           common.legend = TRUE, 
                           legend = "right")

require(grid)

jpeg(file = "Fig1.jpeg", width=960, height=480, units = "px")
annotate_figure(
  fig1, 
  bottom = textGrob("Month"),
  left = textGrob(
    expression(log[10]~(mass~(g))), 
    rot = 90,          # Rotate the y-axis label 90 degrees
    gp = gpar(fontsize = 12)  # Optionally adjust the font size
  )
)

setwd(figures_wd)
jpeg(file = "Fig1.jpeg", width=960, height=480, units = "px")
annotate_figure(
  fig1, 
  bottom = textGrob("Month"),
  left = textGrob(
    expression(log[10]~(mass~(g))), 
    rot = 90,          # Rotate the y-axis label 90 degrees
    gp = gpar(fontsize = 12)  # Optionally adjust the font size
  )
)
dev.off()
## png 
##   2

SMR MODEL

Visualize SMR data

hist(ds1$SMR)        # plot all raw SMR values in a histogram, note skewness

hist(log10(ds1$SMR)) # plot all log-SMR values in a histogram

ggplot(ds1, aes(x = mass, y = SMR, group = ID_fish)) +          # plot raw SMR data by ID against mass
  geom_point(size = 2) +                                        # size of data points
  facet_wrap(. ~ Temp_class) +                                  # produces one plot per temperature treatment
  theme_bw() 

Visualize log-mass and SMR relationship

ggplot(ds1, aes(x = log_mass, y = SMR, group = ID_fish)) +          # plot raw SMR data by ID against log10(mass)
  geom_point(size = 2) +                                            # size of data points
  facet_wrap(. ~ Temp_class) +                                      # produces one plot per temperature treatment
  theme_bw() 

Note that this relationship is non-linear, we expect this to be described by a power function y = a*x^b

This reflects our understanding of the relationship between mass and metabolism, which can be described by the power function y = a*x^b

ggplot(ds1, aes(x = log_mass, y = SMR, group = ID_fish)) +          # plot raw SMR data by ID against log10(mass)
  geom_point(size = 2) +                                            # size of data points
  facet_wrap(. ~ Temp_class) +                                      # produces one plot per temperature treatment
  theme_bw() +
  scale_y_log10()                                                   # plot data with log10 scaled y axis

The trends become linear when SMR is presented on log scale

Univariate SMR models for two temperature treatments

# model structure using log10 SMR and mass for fish at 18C only
SMR18_mod<-lmer(log10(SMR) ~ log_mass + (1 |ID_fish), # random slopes by month, random intercepts by ID_fish
                 data=ds18, REML = T)

summary(SMR18_mod)                                               # generate model output
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(SMR) ~ log_mass + (1 | ID_fish)
##    Data: ds18
## 
## REML criterion at convergence: -553.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6860 -0.7401  0.0462  0.7148  2.7433 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. 
##  ID_fish  (Intercept) 1.063e-17 3.261e-09
##  Residual             9.832e-03 9.916e-02
## Number of obs: 317, groups:  ID_fish, 72
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  -0.67878    0.01415 315.00000  -47.96   <2e-16 ***
## log_mass      0.78058    0.02562 315.00000   30.46   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## log_mass -0.919
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
plot(SMR18_mod)                                                  # plot model residuals

ranova(SMR18_mod)                                                # anova-like test shows insignificant random slopes
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log10(SMR) ~ log_mass + (1 | ID_fish)
##               npar logLik     AIC LRT Df Pr(>Chisq)
## <none>           4 276.78 -545.56                  
## (1 | ID_fish)    3 276.78 -547.56   0  1          1
SMR18_mod_1<-lmer(log10(SMR) ~ log_mass + (1 |ID_fish),          # random slopes term removed
                 data=ds18, REML = T)

summary(SMR18_mod_1)                                             # generate model output
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(SMR) ~ log_mass + (1 | ID_fish)
##    Data: ds18
## 
## REML criterion at convergence: -553.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6860 -0.7401  0.0462  0.7148  2.7433 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. 
##  ID_fish  (Intercept) 1.063e-17 3.261e-09
##  Residual             9.832e-03 9.916e-02
## Number of obs: 317, groups:  ID_fish, 72
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  -0.67878    0.01415 315.00000  -47.96   <2e-16 ***
## log_mass      0.78058    0.02562 315.00000   30.46   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## log_mass -0.919
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
plot(SMR18_mod_1)                                                # plot model residuals

ranova(SMR18_mod_1)                                              # anova-like test shows insignificant random intercepts, but retain in model for testing
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log10(SMR) ~ log_mass + (1 | ID_fish)
##               npar logLik     AIC LRT Df Pr(>Chisq)
## <none>           4 276.78 -545.56                  
## (1 | ID_fish)    3 276.78 -547.56   0  1          1
# model structure using log10 SMR and mass for fish at 23C only
SMR23_mod<-lmer(log10(SMR) ~ log_mass + (1 + Month_int|ID_fish), # random slopes by month, random intercepts by ID_fish
                 data=ds23, REML=T)

summary(SMR23_mod)                                               # generate model output
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(SMR) ~ log_mass + (1 + Month_int | ID_fish)
##    Data: ds23
## 
## REML criterion at convergence: -649.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1396 -0.5503 -0.0023  0.5686  3.2950 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  ID_fish  (Intercept) 0.0052305 0.07232       
##           Month_int   0.0003171 0.01781  -0.94
##  Residual             0.0056041 0.07486       
## Number of obs: 307, groups:  ID_fish, 64
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  -0.47671    0.01133 101.03860  -42.07   <2e-16 ***
## log_mass      0.61389    0.02319 128.47109   26.48   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## log_mass -0.876
plot(SMR23_mod)                                                  # plot model residuals

ranova(SMR23_mod)                                                # anova-like test shows insignificant random slopes
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log10(SMR) ~ log_mass + (1 + Month_int | ID_fish)
##                                        npar logLik     AIC    LRT Df Pr(>Chisq)
## <none>                                    6 324.70 -637.41                     
## Month_int in (1 + Month_int | ID_fish)    4 320.96 -633.92 7.4863  2    0.02368
##                                         
## <none>                                  
## Month_int in (1 + Month_int | ID_fish) *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
SMR23_mod_1<-lmer(log10(SMR) ~ log_mass + (1 |ID_fish),          # random slopes term removed
                 data=ds23, REML=T)

summary(SMR23_mod_1)                                             # generate model output
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(SMR) ~ log_mass + (1 | ID_fish)
##    Data: ds23
## 
## REML criterion at convergence: -641.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0810 -0.5308 -0.0445  0.5962  3.1721 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  ID_fish  (Intercept) 0.0006479 0.02545 
##  Residual             0.0064034 0.08002 
## Number of obs: 307, groups:  ID_fish, 64
## 
## Fixed effects:
##              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)  -0.47203    0.01064 133.59553  -44.38   <2e-16 ***
## log_mass      0.60664    0.02244 178.85157   27.03   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##          (Intr)
## log_mass -0.851
plot(SMR23_mod_1)                                                # plot model residuals

ranova(SMR23_mod_1)                                              # anova-like test shows insignificant random intercepts, but retain in model for testing
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log10(SMR) ~ log_mass + (1 | ID_fish)
##               npar logLik     AIC    LRT Df Pr(>Chisq)  
## <none>           4 320.96 -633.92                       
## (1 | ID_fish)    3 318.78 -631.56 4.3593  1    0.03681 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Combined SMR model

# model using log10 SMR and the interaction of log10 mass and temperature treatment
SMR_mod <- lmer(log10(SMR) ~ log_mass*as.factor(Temp_class) + Sex + Sex*as.factor(Temp_class) + Sex*as.factor(Temp_class)*log_mass + Density + (1|ID_fish),  
                 data=ds1, REML=T)
plot(SMR_mod)                                             # generate model output

summary(SMR_mod)                                          # plot model residuals
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## log10(SMR) ~ log_mass * as.factor(Temp_class) + Sex + Sex * as.factor(Temp_class) +  
##     Sex * as.factor(Temp_class) * log_mass + Density + (1 | ID_fish)
##    Data: ds1
## 
## REML criterion at convergence: -832.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.87616 -0.67366  0.01528  0.66654  2.90845 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  ID_fish  (Intercept) 0.0001859 0.01364 
##  Residual             0.0082037 0.09057 
## Number of obs: 453, groups:  ID_fish, 91
## 
## Fixed effects:
##                                         Estimate Std. Error         df t value
## (Intercept)                            -0.676339   0.031184 156.099873 -21.689
## log_mass                                0.796149   0.032692 255.489558  24.353
## as.factor(Temp_class)23                 0.165816   0.033347 188.402172   4.972
## SexM                                   -0.057479   0.032974 278.565445  -1.743
## Density                                -0.002566   0.003298 152.321130  -0.778
## log_mass:as.factor(Temp_class)23       -0.089839   0.063547 212.044147  -1.414
## as.factor(Temp_class)23:SexM            0.086503   0.047344 220.930300   1.827
## log_mass:SexM                           0.160075   0.063848 346.693411   2.507
## log_mass:as.factor(Temp_class)23:SexM  -0.200908   0.093156 271.605997  -2.157
##                                       Pr(>|t|)    
## (Intercept)                            < 2e-16 ***
## log_mass                               < 2e-16 ***
## as.factor(Temp_class)23               1.48e-06 ***
## SexM                                    0.0824 .  
## Density                                 0.4377    
## log_mass:as.factor(Temp_class)23        0.1589    
## as.factor(Temp_class)23:SexM            0.0690 .  
## log_mass:SexM                           0.0126 *  
## log_mass:as.factor(Temp_class)23:SexM   0.0319 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lg_mss as.(T_)23 SexM   Densty lg_:.(T_)23 a.(T_)23: lg_:SM
## log_mass    -0.663                                                            
## as.fc(T_)23 -0.476  0.565                                                     
## SexM        -0.372  0.559  0.359                                              
## Density     -0.769  0.091  0.122    -0.018                                    
## lg_:.(T_)23  0.390 -0.520 -0.925    -0.286 -0.111                             
## a.(T_)23:SM  0.294 -0.393 -0.698    -0.696 -0.032  0.646                      
## log_mss:SxM  0.290 -0.506 -0.282    -0.925  0.018  0.259       0.644          
## l_:.(T_)23: -0.231  0.351  0.626     0.633  0.030 -0.677      -0.916    -0.685
ranova(SMR_mod)                                           # anova-like test for significance of model components
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log10(SMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 | ID_fish) + log_mass:as.factor(Temp_class) + as.factor(Temp_class):Sex + log_mass:Sex + log_mass:as.factor(Temp_class):Sex
##               npar logLik     AIC     LRT Df Pr(>Chisq)
## <none>          11 416.12 -810.24                      
## (1 | ID_fish)   10 415.93 -811.85 0.38676  1      0.534
# model using log10 SMR and the interaction of log10 mass and temperature treatment with test of random slopes by month
SMR_mod_month <- lmer(log10(SMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (1 + Month_int|ID_fish),  
                 data=ds1, REML=T)
plot(SMR_mod_month)                                             # generate model output

summary(SMR_mod_month)                                          # plot model residuals
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(SMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +  
##     (1 + Month_int | ID_fish)
##    Data: ds1
## 
## REML criterion at convergence: -839.6
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.95365 -0.62500  0.02573  0.65621  2.95222 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  ID_fish  (Intercept) 1.150e-03 0.033906      
##           Month_int   4.241e-05 0.006512 -1.00
##  Residual             8.165e-03 0.090361      
## Number of obs: 453, groups:  ID_fish, 91
## 
## Fixed effects:
##                                    Estimate Std. Error         df t value
## (Intercept)                       -0.693448   0.029255 209.328902 -23.704
## log_mass                           0.835361   0.027904 309.422980  29.937
## as.factor(Temp_class)23            0.195188   0.023171 174.966869   8.424
## SexM                               0.014368   0.009421 121.002750   1.525
## Density                           -0.003163   0.003287 220.153668  -0.962
## log_mass:as.factor(Temp_class)23  -0.154606   0.043270 289.652921  -3.573
##                                  Pr(>|t|)    
## (Intercept)                       < 2e-16 ***
## log_mass                          < 2e-16 ***
## as.factor(Temp_class)23          1.29e-14 ***
## SexM                             0.129854    
## Density                          0.336940    
## log_mass:as.factor(Temp_class)23 0.000413 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lg_mss a.(T_) SexM   Densty
## log_mass    -0.619                            
## as.fc(T_)23 -0.481  0.639                     
## SexM        -0.212  0.215  0.100              
## Density     -0.816  0.113  0.125 -0.020       
## lg_:.(T_)23  0.432 -0.651 -0.915 -0.145 -0.112
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
ranova(SMR_mod_month)                                           # anova-like test for significance of model components
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log10(SMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 + Month_int | ID_fish) + log_mass:as.factor(Temp_class)
##                                        npar logLik     AIC     LRT Df
## <none>                                   10 419.78 -819.56           
## Month_int in (1 + Month_int | ID_fish)    8 419.39 -822.78 0.78507  2
##                                        Pr(>Chisq)
## <none>                                           
## Month_int in (1 + Month_int | ID_fish)     0.6753
# model with random slopes by temp treatment
SMR_mod_trt <- lmer(log10(SMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (0+Temp_class|ID_fish),  
                 data=ds1, REML=T)

plot(SMR_mod_trt)                                             # generate model output

summary(SMR_mod_trt)                                          # plot model residuals
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(SMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +  
##     (0 + Temp_class | ID_fish)
##    Data: ds1
## 
## REML criterion at convergence: -838.8
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.95727 -0.63762  0.00871  0.65226  2.96919 
## 
## Random effects:
##  Groups   Name         Variance  Std.Dev. Corr
##  ID_fish  Temp_class18 4.911e-05 0.007008     
##           Temp_class23 1.882e-04 0.013718 0.56
##  Residual              8.348e-03 0.091368     
## Number of obs: 453, groups:  ID_fish, 91
## 
## Fixed effects:
##                                    Estimate Std. Error         df t value
## (Intercept)                       -0.696017   0.028900 153.536942 -24.084
## log_mass                           0.834470   0.027427 187.550606  30.425
## as.factor(Temp_class)23            0.195939   0.022670 184.796780   8.643
## SexM                               0.015418   0.009284  87.658721   1.661
## Density                           -0.002715   0.003269 150.116946  -0.831
## log_mass:as.factor(Temp_class)23  -0.154162   0.042967 211.007690  -3.588
##                                  Pr(>|t|)    
## (Intercept)                       < 2e-16 ***
## log_mass                          < 2e-16 ***
## as.factor(Temp_class)23          2.56e-15 ***
## SexM                             0.100349    
## Density                          0.407416    
## log_mass:as.factor(Temp_class)23 0.000414 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lg_mss a.(T_) SexM   Densty
## log_mass    -0.609                            
## as.fc(T_)23 -0.471  0.632                     
## SexM        -0.211  0.218  0.096              
## Density     -0.823  0.112  0.126 -0.020       
## lg_:.(T_)23  0.420 -0.643 -0.912 -0.141 -0.108
ranova(SMR_mod_trt)                                           # anova-like test for significance of model components
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log10(SMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (0 + Temp_class | ID_fish) + log_mass:as.factor(Temp_class)
##                                          npar logLik     AIC      LRT Df
## <none>                                     10 419.42 -818.84            
## Temp_class in (0 + Temp_class | ID_fish)    8 419.39 -822.78 0.061499  2
##                                          Pr(>Chisq)
## <none>                                             
## Temp_class in (0 + Temp_class | ID_fish)     0.9697
# model with random intercepts by temp treatment
SMR_mod_trtint <- lmer(log10(SMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (1|Temp_class),  
                 data=ds1, REML=T)

plot(SMR_mod_trtint)                                             # generate model output

summary(SMR_mod_trtint)                                          # plot model residuals
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(SMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +  
##     (1 | Temp_class)
##    Data: ds1
## 
## REML criterion at convergence: -838.7
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.96092 -0.63658  0.00346  0.65601  2.98231 
## 
## Random effects:
##  Groups     Name        Variance  Std.Dev.
##  Temp_class (Intercept) 0.0002785 0.01669 
##  Residual               0.0084475 0.09191 
## Number of obs: 453, groups:  Temp_class, 2
## 
## Fixed effects:
##                                    Estimate Std. Error         df t value
## (Intercept)                      -6.951e-01  3.302e-02  2.625e-08 -21.051
## log_mass                          8.339e-01  2.736e-02  4.470e+02  30.485
## as.factor(Temp_class)23           1.941e-01  3.247e-02  6.139e-09   5.978
## SexM                              1.524e-02  9.057e-03  4.470e+02   1.682
## Density                          -2.806e-03  3.211e-03  4.470e+02  -0.874
## log_mass:as.factor(Temp_class)23 -1.502e-01  4.233e-02  4.470e+02  -3.548
##                                  Pr(>|t|)    
## (Intercept)                      1.000000    
## log_mass                          < 2e-16 ***
## as.factor(Temp_class)23          1.000000    
## SexM                             0.093181 .  
## Density                          0.382747    
## log_mass:as.factor(Temp_class)23 0.000429 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lg_mss a.(T_) SexM   Densty
## log_mass    -0.528                            
## as.fc(T_)23 -0.545  0.440                     
## SexM        -0.180  0.216  0.067              
## Density     -0.707  0.109  0.088 -0.021       
## lg_:.(T_)23  0.371 -0.651 -0.629 -0.141 -0.111
ranova(SMR_mod_trtint)                                           # anova-like test for significance of model components
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log10(SMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 | Temp_class) + log_mass:as.factor(Temp_class)
##                  npar logLik     AIC         LRT Df Pr(>Chisq)
## <none>              8 419.33 -822.66                          
## (1 | Temp_class)    7 419.33 -824.66 -1.1369e-13  1          1
# test sex terms for model inclusion

SMR_mod_sex <- lmer(log10(SMR) ~ 0 + as.factor(Sex):as.factor(Temp_class) + as.factor(Sex):as.factor(Temp_class):log_mass + Density + (1 |ID_fish),  
                 data=ds1, REML=T)
plot(SMR_mod_sex)

summary(SMR_mod_sex)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## log10(SMR) ~ 0 + as.factor(Sex):as.factor(Temp_class) + as.factor(Sex):as.factor(Temp_class):log_mass +  
##     Density + (1 | ID_fish)
##    Data: ds1
## 
## REML criterion at convergence: -832.2
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.87616 -0.67366  0.01528  0.66654  2.90845 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  ID_fish  (Intercept) 0.0001859 0.01364 
##  Residual             0.0082037 0.09057 
## Number of obs: 453, groups:  ID_fish, 91
## 
## Fixed effects:
##                                                    Estimate Std. Error
## Density                                           -0.002566   0.003298
## as.factor(Sex)F:as.factor(Temp_class)18           -0.676339   0.031184
## as.factor(Sex)M:as.factor(Temp_class)18           -0.733818   0.035973
## as.factor(Sex)F:as.factor(Temp_class)23           -0.510523   0.033096
## as.factor(Sex)M:as.factor(Temp_class)23           -0.481499   0.030654
## as.factor(Sex)F:as.factor(Temp_class)18:log_mass   0.796149   0.032692
## as.factor(Sex)M:as.factor(Temp_class)18:log_mass   0.956224   0.055066
## as.factor(Sex)F:as.factor(Temp_class)23:log_mass   0.706309   0.054270
## as.factor(Sex)M:as.factor(Temp_class)23:log_mass   0.665476   0.040820
##                                                          df t value Pr(>|t|)
## Density                                          152.321130  -0.778    0.438
## as.factor(Sex)F:as.factor(Temp_class)18          156.099871 -21.689   <2e-16
## as.factor(Sex)M:as.factor(Temp_class)18          226.413501 -20.399   <2e-16
## as.factor(Sex)F:as.factor(Temp_class)23          149.558833 -15.426   <2e-16
## as.factor(Sex)M:as.factor(Temp_class)23          158.947247 -15.707   <2e-16
## as.factor(Sex)F:as.factor(Temp_class)18:log_mass 255.489556  24.353   <2e-16
## as.factor(Sex)M:as.factor(Temp_class)18:log_mass 374.137705  17.365   <2e-16
## as.factor(Sex)F:as.factor(Temp_class)23:log_mass 198.971185  13.015   <2e-16
## as.factor(Sex)M:as.factor(Temp_class)23:log_mass 256.884341  16.303   <2e-16
##                                                     
## Density                                             
## as.factor(Sex)F:as.factor(Temp_class)18          ***
## as.factor(Sex)M:as.factor(Temp_class)18          ***
## as.factor(Sex)F:as.factor(Temp_class)23          ***
## as.factor(Sex)M:as.factor(Temp_class)23          ***
## as.factor(Sex)F:as.factor(Temp_class)18:log_mass ***
## as.factor(Sex)M:as.factor(Temp_class)18:log_mass ***
## as.factor(Sex)F:as.factor(Temp_class)23:log_mass ***
## as.factor(Sex)M:as.factor(Temp_class)23:log_mass ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##                 Densty as.(S)F:.(T_)18 as.(S)M:.(T_)18 as.(S)F:.(T_)23
## as.(S)F:.(T_)18 -0.769                                                
## as.(S)M:.(T_)18 -0.683  0.526                                         
## as.(S)F:.(T_)23 -0.602  0.463           0.411                         
## as.(S)M:.(T_)23 -0.719  0.553           0.491           0.433         
## a.(S)F:.(T_)18:  0.091 -0.663          -0.062          -0.055         
## a.(S)M:.(T_)18:  0.074 -0.057          -0.729          -0.045         
## a.(S)F:.(T_)23: -0.075  0.058           0.051          -0.694         
## a.(S)M:.(T_)23: -0.004  0.003           0.003           0.002         
##                 as.(S)M:.(T_)23 a.(S)F:.(T_)18: a.(S)M:.(T_)18: a.(S)F:.(T_)23:
## as.(S)F:.(T_)18                                                                
## as.(S)M:.(T_)18                                                                
## as.(S)F:.(T_)23                                                                
## as.(S)M:.(T_)23                                                                
## a.(S)F:.(T_)18: -0.065                                                         
## a.(S)M:.(T_)18: -0.054           0.007                                         
## a.(S)F:.(T_)23:  0.054          -0.007          -0.006                         
## a.(S)M:.(T_)23: -0.604           0.000           0.000           0.000

Again, our ranova analysis shows that the random intercepts by ID_fish are not significant to this model. However, this will be retained because repeated measures were made on the same individuals, and so that the model can be included in the below covariance matrix.

Plot SMR vs. log mass mean predicted level

pred_SMR <- ggpredict(SMR_mod, terms = c("log_mass", "Temp_class"))

ggplot() +
  geom_errorbar(data = pred_SMR, aes(x=x, ymin = conf.low, ymax = conf.high, color = group, group = group), width = 0.2) +   # add error bars showing 95% confidence intervals
  geom_point(data = pred_SMR, aes(x = x, y = predicted, color = group, group = group)) +  
  geom_line(data = pred_SMR, aes(x = x, y = predicted, color = group, group = group)) +
  labs(x = "log_mass", y = "Predicted SMR", color = "Temp_class") +
  theme_bw() +
  scale_color_manual(values = c("#78B7C5", "#F21A00")) +
  scale_y_log10() 

F test on interaction effect

anova(SMR_mod, type = "3")
## Type III Analysis of Variance Table with Satterthwaite's method
##                                    Sum Sq Mean Sq NumDF  DenDF   F value
## log_mass                           9.2265  9.2265     1 269.79 1124.6757
## as.factor(Temp_class)              0.6283  0.6283     1 222.17   76.5912
## Sex                                0.0030  0.0030     1 223.64    0.3604
## Density                            0.0050  0.0050     1 152.32    0.6056
## log_mass:as.factor(Temp_class)     0.1351  0.1351     1 273.42   16.4634
## as.factor(Temp_class):Sex          0.0274  0.0274     1 220.93    3.3384
## log_mass:Sex                       0.0134  0.0134     1 274.33    1.6352
## log_mass:as.factor(Temp_class):Sex 0.0382  0.0382     1 271.61    4.6513
##                                       Pr(>F)    
## log_mass                           < 2.2e-16 ***
## as.factor(Temp_class)              5.351e-16 ***
## Sex                                  0.54888    
## Density                              0.43767    
## log_mass:as.factor(Temp_class)     6.477e-05 ***
## as.factor(Temp_class):Sex            0.06903 .  
## log_mass:Sex                         0.20207    
## log_mass:as.factor(Temp_class):Sex   0.03191 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Generate fitted data

z_SMR<-augment(SMR_mod)                             # get variables and fitted values from model
head(z_SMR)
## # A tibble: 6 × 18
##   .rownames `log10(SMR)` log_mass `as.factor(Temp_class)` Sex   Density ID_fish
##   <chr>            <dbl>    <dbl> <fct>                   <fct>   <dbl>   <dbl>
## 1 1              -0.0733   0.423  23                      F           8      81
## 2 2              -0.546   -0.0915 23                      M           8      82
## 3 3              -0.276    0.250  23                      F           8      83
## 4 4              -0.0948   0.436  23                      M           8      84
## 5 5              -0.151    0.360  23                      M           8      85
## 6 7              -0.296    0.310  23                      M           8      87
## # ℹ 11 more variables: .fitted <dbl>, .resid <dbl>, .hat <dbl>, .cooksd <dbl>,
## #   .fixed <dbl>, .mu <dbl>, .offset <dbl>, .sqrtXwt <dbl>, .sqrtrwt <dbl>,
## #   .weights <dbl>, .wtres <dbl>
z_SMR <- z_SMR %>% 
  mutate(
    pwr.fitted = 10^(.fitted)
  )
##### check assumptions: plot all residual values against the fitted 
ggplot(z_SMR, aes(x = .fitted, y = .resid)) +   # same plot as before
  geom_point(size = 2) 

ggplot(z_SMR, aes(x = log_mass, y = SMR, group = ID_fish)) + 
  geom_line(aes(y = pwr.fitted, color = `as.factor(Temp_class)`),      # note this will look odd if there are additional fixed effects      
            alpha = 0.6) +
  theme_classic() +
  scale_color_manual(values  = c("#78B7C5", "#F21A00")) +
  scale_y_log10()+
  labs(x = expression(log[10]~(mass~(g))), y = expression(Fitted~SMR~(mg~O[2]~h^-1)), color = expression("Temperature ("*degree*C*")"))

At this stage, we have found a strong model for SMR which captures variance by treatment and uses centered mass. The next step is to repeat this modeling approach for RMR and MMR.

RMR MODEL

# model using log10 RMR and the interaction of log10 mass and temperature treatment
RMR_mod <- lmer(log10(RMR) ~ log_mass*as.factor(Temp_class) + Sex + Sex*as.factor(Temp_class) + Sex*as.factor(Temp_class)*log_mass + Density + (1 |ID_fish),  
                 data=ds1, REML=T)
plot(RMR_mod)

summary(RMR_mod)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## log10(RMR) ~ log_mass * as.factor(Temp_class) + Sex + Sex * as.factor(Temp_class) +  
##     Sex * as.factor(Temp_class) * log_mass + Density + (1 | ID_fish)
##    Data: ds1
## 
## REML criterion at convergence: -895.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6217 -0.6035  0.0362  0.6597  2.4226 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID_fish  (Intercept) 0.001370 0.03701 
##  Residual             0.006309 0.07943 
## Number of obs: 453, groups:  ID_fish, 91
## 
## Fixed effects:
##                                         Estimate Std. Error         df t value
## (Intercept)                            -0.458863   0.033569 211.189136 -13.669
## log_mass                                0.715302   0.032984 359.999351  21.686
## as.factor(Temp_class)23                 0.133553   0.035094 249.693441   3.806
## SexM                                   -0.053158   0.033240 333.370196  -1.599
## Density                                -0.004299   0.003555 213.356425  -1.209
## log_mass:as.factor(Temp_class)23       -0.097402   0.065508 312.941087  -1.487
## as.factor(Temp_class)23:SexM            0.070119   0.049047 273.943840   1.430
## log_mass:SexM                           0.124388   0.062014 428.861974   2.006
## log_mass:as.factor(Temp_class)23:SexM  -0.163282   0.093374 373.092786  -1.749
##                                       Pr(>|t|)    
## (Intercept)                            < 2e-16 ***
## log_mass                               < 2e-16 ***
## as.factor(Temp_class)23               0.000178 ***
## SexM                                  0.110717    
## Density                               0.227872    
## log_mass:as.factor(Temp_class)23      0.138054    
## as.factor(Temp_class)23:SexM          0.153960    
## log_mass:SexM                         0.045505 *  
## log_mass:as.factor(Temp_class)23:SexM 0.081167 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lg_mss as.(T_)23 SexM   Densty lg_:.(T_)23 a.(T_)23: lg_:SM
## log_mass    -0.655                                                            
## as.fc(T_)23 -0.456  0.542                                                     
## SexM        -0.381  0.555  0.362                                              
## Density     -0.791  0.134  0.124     0.004                                    
## lg_:.(T_)23  0.363 -0.509 -0.904    -0.280 -0.110                             
## a.(T_)23:SM  0.289 -0.382 -0.710    -0.678 -0.042  0.642                      
## log_mss:SxM  0.295 -0.523 -0.280    -0.895 -0.003  0.263       0.607          
## l_:.(T_)23: -0.224  0.352  0.629     0.595  0.038 -0.697      -0.887    -0.664
ranova(RMR_mod)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log10(RMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 | ID_fish) + log_mass:as.factor(Temp_class) + as.factor(Temp_class):Sex + log_mass:Sex + log_mass:as.factor(Temp_class):Sex
##               npar logLik     AIC    LRT Df Pr(>Chisq)    
## <none>          11 447.95 -873.91                         
## (1 | ID_fish)   10 437.96 -855.92 19.988  1  7.792e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# model using log10 SMR and the interaction of log10 mass and temperature treatment with test of random slopes by month
RMR_mod_month <- lmer(log10(RMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (1 + Month_int|ID_fish),  
                 data=ds1, REML=T)
plot(RMR_mod_month)                                             # generate model output

summary(RMR_mod_month)                                          # plot model residuals
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(RMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +  
##     (1 + Month_int | ID_fish)
##    Data: ds1
## 
## REML criterion at convergence: -904.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6681 -0.5983  0.0284  0.6278  2.4239 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  ID_fish  (Intercept) 1.509e-03 0.038848      
##           Month_int   1.135e-06 0.001065 -1.00
##  Residual             6.380e-03 0.079873      
## Number of obs: 453, groups:  ID_fish, 91
## 
## Fixed effects:
##                                    Estimate Std. Error         df t value
## (Intercept)                       -0.476881   0.031184 208.946377 -15.293
## log_mass                           0.748619   0.027720 371.916742  27.006
## as.factor(Temp_class)23            0.159088   0.023714 198.905581   6.709
## SexM                               0.003107   0.011023  86.886610   0.282
## Density                           -0.004258   0.003515 203.921547  -1.211
## log_mass:as.factor(Temp_class)23  -0.153474   0.043242 358.751139  -3.549
##                                  Pr(>|t|)    
## (Intercept)                       < 2e-16 ***
## log_mass                          < 2e-16 ***
## as.factor(Temp_class)23             2e-10 ***
## SexM                             0.778752    
## Density                          0.227163    
## log_mass:as.factor(Temp_class)23 0.000438 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lg_mss a.(T_) SexM   Densty
## log_mass    -0.610                            
## as.fc(T_)23 -0.458  0.615                     
## SexM        -0.203  0.178  0.064              
## Density     -0.836  0.152  0.128 -0.012       
## lg_:.(T_)23  0.404 -0.644 -0.886 -0.118 -0.112
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 0.00514734 (tol = 0.002, component 1)
ranova(RMR_mod_month)                                           # anova-like test for significance of model components
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log10(RMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 + Month_int | ID_fish) + log_mass:as.factor(Temp_class)
##                                        npar logLik     AIC      LRT Df
## <none>                                   10 452.28 -884.56            
## Month_int in (1 + Month_int | ID_fish)    8 452.24 -888.49 0.076141  2
##                                        Pr(>Chisq)
## <none>                                           
## Month_int in (1 + Month_int | ID_fish)     0.9626
# model with random intercepts by temp treatment
RMR_mod_trtint <- lmer(log10(RMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (1|Temp_class),  
                 data=ds1, REML=T)
plot(RMR_mod_trtint)

summary(RMR_mod_trtint)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(RMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +  
##     (1 | Temp_class)
##    Data: ds1
## 
## REML criterion at convergence: -886.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0205 -0.6163  0.0171  0.6597  3.1772 
## 
## Random effects:
##  Groups     Name        Variance  Std.Dev.
##  Temp_class (Intercept) 9.358e-05 0.009674
##  Residual               7.588e-03 0.087111
## Number of obs: 453, groups:  Temp_class, 2
## 
## Fixed effects:
##                                    Estimate Std. Error         df t value
## (Intercept)                      -4.775e-01  2.868e-02  3.908e-08 -16.649
## log_mass                          7.495e-01  2.593e-02  4.470e+02  28.908
## as.factor(Temp_class)23           1.463e-01  2.518e-02  5.802e-09   5.811
## SexM                              3.679e-03  8.584e-03  4.470e+02   0.429
## Density                          -4.302e-03  3.044e-03  4.470e+02  -1.413
## log_mass:as.factor(Temp_class)23 -1.232e-01  4.012e-02  4.470e+02  -3.071
##                                  Pr(>|t|)    
## (Intercept)                       1.00000    
## log_mass                          < 2e-16 ***
## as.factor(Temp_class)23           1.00000    
## SexM                              0.66846    
## Density                           0.15824    
## log_mass:as.factor(Temp_class)23  0.00226 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lg_mss a.(T_) SexM   Densty
## log_mass    -0.576                            
## as.fc(T_)23 -0.510  0.538                     
## SexM        -0.197  0.216  0.081              
## Density     -0.772  0.109  0.108 -0.021       
## lg_:.(T_)23  0.404 -0.651 -0.768 -0.141 -0.111
ranova(RMR_mod_trtint)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log10(RMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 | Temp_class) + log_mass:as.factor(Temp_class)
##                  npar logLik    AIC        LRT Df Pr(>Chisq)
## <none>              8  443.3 -870.6                         
## (1 | Temp_class)    7  443.3 -872.6 2.2737e-13  1          1
# model with random slopes by temp treatment
RMR_mod_trt <- lmer(log10(RMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (0 + Temp_class|ID_fish),  
                 data=ds1, REML=T)
plot(RMR_mod_trt)

summary(RMR_mod_trt)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(RMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +  
##     (0 + Temp_class | ID_fish)
##    Data: ds1
## 
## REML criterion at convergence: -904.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6323 -0.5883  0.0271  0.6401  2.3369 
## 
## Random effects:
##  Groups   Name         Variance Std.Dev. Corr
##  ID_fish  Temp_class18 0.001392 0.03731      
##           Temp_class23 0.001059 0.03254  0.15
##  Residual              0.006388 0.07992      
## Number of obs: 453, groups:  ID_fish, 91
## 
## Fixed effects:
##                                    Estimate Std. Error         df t value
## (Intercept)                       -0.476861   0.031190 208.937568 -15.289
## log_mass                           0.748448   0.027931 302.207977  26.796
## as.factor(Temp_class)23            0.158264   0.023388 232.444603   6.767
## SexM                               0.003190   0.010952  82.507699   0.291
## Density                           -0.004258   0.003507 196.713602  -1.214
## log_mass:as.factor(Temp_class)23  -0.150036   0.042951 298.007754  -3.493
##                                  Pr(>|t|)    
## (Intercept)                       < 2e-16 ***
## log_mass                          < 2e-16 ***
## as.factor(Temp_class)23          1.06e-10 ***
## SexM                              0.77155    
## Density                           0.22619    
## log_mass:as.factor(Temp_class)23  0.00055 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lg_mss a.(T_) SexM   Densty
## log_mass    -0.607                            
## as.fc(T_)23 -0.470  0.624                     
## SexM        -0.197  0.169  0.057              
## Density     -0.835  0.151  0.135 -0.013       
## lg_:.(T_)23  0.410 -0.653 -0.886 -0.112 -0.117
## optimizer (nloptwrap) convergence code: 0 (OK)
## unable to evaluate scaled gradient
## Model failed to converge: degenerate  Hessian with 1 negative eigenvalues
ranova(RMR_mod_trt)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log10(RMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (0 + Temp_class | ID_fish) + log_mass:as.factor(Temp_class)
##                                          npar logLik     AIC     LRT Df
## <none>                                     10 452.33 -884.66           
## Temp_class in (0 + Temp_class | ID_fish)    8 452.24 -888.49 0.17085  2
##                                          Pr(>Chisq)
## <none>                                             
## Temp_class in (0 + Temp_class | ID_fish)     0.9181
# test sex terms for model inclusion

RMR_mod_sex <- lmer(log10(RMR) ~ 0 + as.factor(Sex):as.factor(Temp_class) + as.factor(Sex):as.factor(Temp_class):log_mass + Density + (1 |ID_fish),  
                 data=ds1, REML=T)
plot(RMR_mod_sex)

summary(RMR_mod_sex)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## log10(RMR) ~ 0 + as.factor(Sex):as.factor(Temp_class) + as.factor(Sex):as.factor(Temp_class):log_mass +  
##     Density + (1 | ID_fish)
##    Data: ds1
## 
## REML criterion at convergence: -895.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.6217 -0.6035  0.0362  0.6597  2.4226 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  ID_fish  (Intercept) 0.001370 0.03701 
##  Residual             0.006309 0.07943 
## Number of obs: 453, groups:  ID_fish, 91
## 
## Fixed effects:
##                                                    Estimate Std. Error
## Density                                           -0.004299   0.003555
## as.factor(Sex)F:as.factor(Temp_class)18           -0.458863   0.033569
## as.factor(Sex)M:as.factor(Temp_class)18           -0.512022   0.037156
## as.factor(Sex)F:as.factor(Temp_class)23           -0.325311   0.035831
## as.factor(Sex)M:as.factor(Temp_class)23           -0.308350   0.032984
## as.factor(Sex)F:as.factor(Temp_class)18:log_mass   0.715302   0.032984
## as.factor(Sex)M:as.factor(Temp_class)18:log_mass   0.839690   0.052870
## as.factor(Sex)F:as.factor(Temp_class)23:log_mass   0.617900   0.056382
## as.factor(Sex)M:as.factor(Temp_class)23:log_mass   0.579006   0.041101
##                                                          df t value Pr(>|t|)
## Density                                          213.356426  -1.209    0.228
## as.factor(Sex)F:as.factor(Temp_class)18          211.189136 -13.669   <2e-16
## as.factor(Sex)M:as.factor(Temp_class)18          283.322955 -13.780   <2e-16
## as.factor(Sex)F:as.factor(Temp_class)23          200.632088  -9.079   <2e-16
## as.factor(Sex)M:as.factor(Temp_class)23          207.520955  -9.348   <2e-16
## as.factor(Sex)F:as.factor(Temp_class)18:log_mass 359.999351  21.686   <2e-16
## as.factor(Sex)M:as.factor(Temp_class)18:log_mass 439.883035  15.882   <2e-16
## as.factor(Sex)F:as.factor(Temp_class)23:log_mass 294.895487  10.959   <2e-16
## as.factor(Sex)M:as.factor(Temp_class)23:log_mass 367.788434  14.088   <2e-16
##                                                     
## Density                                             
## as.factor(Sex)F:as.factor(Temp_class)18          ***
## as.factor(Sex)M:as.factor(Temp_class)18          ***
## as.factor(Sex)F:as.factor(Temp_class)23          ***
## as.factor(Sex)M:as.factor(Temp_class)23          ***
## as.factor(Sex)F:as.factor(Temp_class)18:log_mass ***
## as.factor(Sex)M:as.factor(Temp_class)18:log_mass ***
## as.factor(Sex)F:as.factor(Temp_class)23:log_mass ***
## as.factor(Sex)M:as.factor(Temp_class)23:log_mass ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##                 Densty as.(S)F:.(T_)18 as.(S)M:.(T_)18 as.(S)F:.(T_)23
## as.(S)F:.(T_)18 -0.791                                                
## as.(S)M:.(T_)18 -0.711  0.562                                         
## as.(S)F:.(T_)23 -0.620  0.490           0.441                         
## as.(S)M:.(T_)23 -0.731  0.578           0.520           0.453         
## a.(S)F:.(T_)18:  0.134 -0.655          -0.095          -0.083         
## a.(S)M:.(T_)18:  0.080 -0.063          -0.687          -0.050         
## a.(S)F:.(T_)23: -0.049  0.039           0.035          -0.681         
## a.(S)M:.(T_)23:  0.015 -0.012          -0.011          -0.010         
##                 as.(S)M:.(T_)23 a.(S)F:.(T_)18: a.(S)M:.(T_)18: a.(S)F:.(T_)23:
## as.(S)F:.(T_)18                                                                
## as.(S)M:.(T_)18                                                                
## as.(S)F:.(T_)23                                                                
## as.(S)M:.(T_)23                                                                
## a.(S)F:.(T_)18: -0.098                                                         
## a.(S)M:.(T_)18: -0.059           0.011                                         
## a.(S)F:.(T_)23:  0.036          -0.007          -0.004                         
## a.(S)M:.(T_)23: -0.579           0.002           0.001          -0.001
# plot RMR vs. log mass mean predicted level 

pred_RMR <- ggpredict(RMR_mod, terms = c("log_mass", "Temp_class"))

ggplot() +
  geom_errorbar(data = pred_RMR, aes(x=x, ymin = conf.low, ymax = conf.high, color = group, group = group), width = 0.2) +   # add error bars showing 95% confidence intervals
  geom_point(data = pred_RMR, aes(x = x, y = predicted, color = group, group = group)) +  
  geom_line(data = pred_RMR, aes(x = x, y = predicted, color = group, group = group)) +
  labs(x = "log_mass_centre", y = "Predicted RMR", color = "Temp_class") +
  theme_bw() +
  scale_color_manual(values = c("#78B7C5", "#F21A00")) +
  scale_y_log10() 

# plot the raw data as a cloud around the predictions

ggplot() +
  geom_point(data = ds1, aes(x = log_mass, y = RMR, group = factor(Temp_class), color = factor(Temp_class)), 
             alpha = 0.3, position = position_jitter(width = 0.1)) +
  geom_point(data = pred_RMR, aes(x = x, y = predicted, color = group, group = group), 
             size = 3) +
  geom_line(data = pred_RMR, aes(x = x, y = predicted, color = group, group = group), 
            size = 2) +
  labs(x = "log_mass_centre", y = "Predicted RMR", color = "Temp_class") +
  theme_bw() +
  scale_color_manual(values = c("#78B7C5", "#F21A00")) +
  scale_y_log10() 

# F test on interaction effect

anova(RMR_mod, type = "3")
## Type III Analysis of Variance Table with Satterthwaite's method
##                                    Sum Sq Mean Sq NumDF  DenDF  F value
## log_mass                           5.4611  5.4611     1 367.11 865.6183
## as.factor(Temp_class)              0.2934  0.2934     1 276.53  46.5051
## Sex                                0.0034  0.0034     1 276.43   0.5449
## Density                            0.0092  0.0092     1 213.36   1.4625
## log_mass:as.factor(Temp_class)     0.0917  0.0917     1 375.49  14.5343
## as.factor(Temp_class):Sex          0.0129  0.0129     1 273.94   2.0439
## log_mass:Sex                       0.0053  0.0053     1 375.33   0.8386
## log_mass:as.factor(Temp_class):Sex 0.0193  0.0193     1 373.09   3.0579
##                                       Pr(>F)    
## log_mass                           < 2.2e-16 ***
## as.factor(Temp_class)              5.738e-11 ***
## Sex                                0.4610311    
## Density                            0.2278721    
## log_mass:as.factor(Temp_class)     0.0001608 ***
## as.factor(Temp_class):Sex          0.1539602    
## log_mass:Sex                       0.3603912    
## log_mass:as.factor(Temp_class):Sex 0.0811666 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# generate fitted data

z_RMR<-augment(RMR_mod)                                            # get variables and fitted values from model
head(z_RMR)
## # A tibble: 6 × 18
##   .rownames `log10(RMR)` log_mass `as.factor(Temp_class)` Sex   Density ID_fish
##   <chr>            <dbl>    <dbl> <fct>                   <fct>   <dbl>   <dbl>
## 1 1             -0.0212    0.423  23                      F           8      81
## 2 2             -0.386    -0.0915 23                      M           8      82
## 3 3             -0.147     0.250  23                      F           8      83
## 4 4              0.00674   0.436  23                      M           8      84
## 5 5             -0.0573    0.360  23                      M           8      85
## 6 7             -0.161     0.310  23                      M           8      87
## # ℹ 11 more variables: .fitted <dbl>, .resid <dbl>, .hat <dbl>, .cooksd <dbl>,
## #   .fixed <dbl>, .mu <dbl>, .offset <dbl>, .sqrtXwt <dbl>, .sqrtrwt <dbl>,
## #   .weights <dbl>, .wtres <dbl>
z_RMR <- z_RMR %>% 
  mutate(
    pwr.fitted = 10^(.fitted)
  )

# check assumptions: plot all residual values against the fitted 

ggplot(z_RMR, aes(x = .fitted, y = .resid)) +                      # same plot as before
  geom_point(size = 2) 

ggplot(z_RMR, aes(x = log_mass, y = RMR, group = ID_fish)) + 
  geom_line(aes(y = pwr.fitted, color = `as.factor(Temp_class)`),      
            alpha = 0.6) +
  theme_bw() +
  scale_y_log10() +
  scale_color_manual(values  = c("#78B7C5", "#F21A00"))

MMR MODEL

# model using log10 MMR and the interaction of log10 mass and temperature treatment

MMR_mod <- lmer(log10(MMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (1 |ID_fish),  
                 data=ds1, REML=T)

plot(MMR_mod)

summary(MMR_mod) 
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(MMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +  
##     (1 | ID_fish)
##    Data: ds1
## 
## REML criterion at convergence: -1086.6
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0358 -0.6335  0.0066  0.5706  3.5127 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  ID_fish  (Intercept) 0.0008517 0.02918 
##  Residual             0.0042430 0.06514 
## Number of obs: 453, groups:  ID_fish, 91
## 
## Fixed effects:
##                                    Estimate Std. Error         df t value
## (Intercept)                       -0.034542   0.025478 210.005598  -1.356
## log_mass                           0.718537   0.022624 372.713216  31.760
## as.factor(Temp_class)23            0.050788   0.019226 264.403760   2.642
## SexM                               0.026480   0.009007  86.996070   2.940
## Density                           -0.012575   0.002877 211.708210  -4.371
## log_mass:as.factor(Temp_class)23  -0.096076   0.035291 359.553907  -2.722
##                                  Pr(>|t|)    
## (Intercept)                       0.17663    
## log_mass                          < 2e-16 ***
## as.factor(Temp_class)23           0.00874 ** 
## SexM                              0.00421 ** 
## Density                          1.94e-05 ***
## log_mass:as.factor(Temp_class)23  0.00680 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lg_mss a.(T_) SexM   Densty
## log_mass    -0.605                            
## as.fc(T_)23 -0.457  0.614                     
## SexM        -0.202  0.176  0.060              
## Density     -0.838  0.152  0.130 -0.013       
## lg_:.(T_)23  0.400 -0.644 -0.884 -0.115 -0.111
ranova(MMR_mod)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log10(MMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 | ID_fish) + log_mass:as.factor(Temp_class)
##               npar logLik     AIC    LRT Df Pr(>Chisq)    
## <none>           8 543.31 -1070.6                         
## (1 | ID_fish)    7 534.20 -1054.4 18.232  1  1.956e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# model using log10 SMR and the interaction of log10 mass and temperature treatment with test of random slopes by month
MMR_mod_month <- lmer(log10(MMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (1 + Month_int|ID_fish),  
                 data=ds1, REML=T)
plot(MMR_mod_month)                                             # generate model output

summary(MMR_mod_month)                                          # plot model residuals
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(MMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +  
##     (1 + Month_int | ID_fish)
##    Data: ds1
## 
## REML criterion at convergence: -1090.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.6553 -0.6206  0.0243  0.5554  3.5381 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. Corr 
##  ID_fish  (Intercept) 0.0013142 0.03625       
##           Month_int   0.0001224 0.01106  -0.62
##  Residual             0.0039368 0.06274       
## Number of obs: 453, groups:  ID_fish, 91
## 
## Fixed effects:
##                                    Estimate Std. Error         df t value
## (Intercept)                       -0.037512   0.025652 199.690054  -1.462
## log_mass                           0.716517   0.023329 271.438607  30.714
## as.factor(Temp_class)23            0.051810   0.019205 169.950435   2.698
## SexM                               0.028424   0.008949  87.097271   3.176
## Density                           -0.012195   0.002938 191.627985  -4.151
## log_mass:as.factor(Temp_class)23  -0.096354   0.036605 238.653992  -2.632
##                                  Pr(>|t|)    
## (Intercept)                       0.14523    
## log_mass                          < 2e-16 ***
## as.factor(Temp_class)23           0.00768 ** 
## SexM                              0.00206 ** 
## Density                          4.99e-05 ***
## log_mass:as.factor(Temp_class)23  0.00903 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lg_mss a.(T_) SexM   Densty
## log_mass    -0.576                            
## as.fc(T_)23 -0.449  0.613                     
## SexM        -0.195  0.178  0.055              
## Density     -0.840  0.116  0.123 -0.019       
## lg_:.(T_)23  0.382 -0.639 -0.885 -0.108 -0.093
ranova(MMR_mod_month)                                           # anova-like test for significance of model components
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log10(MMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 + Month_int | ID_fish) + log_mass:as.factor(Temp_class)
##                                        npar logLik     AIC    LRT Df Pr(>Chisq)
## <none>                                   10 545.27 -1070.5                     
## Month_int in (1 + Month_int | ID_fish)    8 543.31 -1070.6 3.9032  2      0.142
# model with random intercepts by temp treatment
MMR_mod_trtint <- lmer(log10(MMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (1|Temp_class),  
                 data=ds1, REML=T)

plot(MMR_mod_trtint)

summary(MMR_mod_trtint) 
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(MMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +  
##     (1 | Temp_class)
##    Data: ds1
## 
## REML criterion at convergence: -1068.4
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.3156 -0.6143  0.0195  0.6269  3.1035 
## 
## Random effects:
##  Groups     Name        Variance  Std.Dev.
##  Temp_class (Intercept) 6.744e-05 0.008212
##  Residual               5.053e-03 0.071082
## Number of obs: 453, groups:  Temp_class, 2
## 
## Fixed effects:
##                                    Estimate Std. Error         df t value
## (Intercept)                      -5.223e-02  2.351e-02  1.327e-07  -2.221
## log_mass                          7.324e-01  2.116e-02  4.470e+02  34.618
## as.factor(Temp_class)23           5.136e-02  2.079e-02  2.029e-08   2.470
## SexM                              2.738e-02  7.005e-03  4.470e+02   3.909
## Density                          -1.111e-02  2.484e-03  4.470e+02  -4.473
## log_mass:as.factor(Temp_class)23 -9.468e-02  3.274e-02  4.470e+02  -2.892
##                                  Pr(>|t|)    
## (Intercept)                      0.999999    
## log_mass                          < 2e-16 ***
## as.factor(Temp_class)23          1.000000    
## SexM                             0.000107 ***
## Density                          9.79e-06 ***
## log_mass:as.factor(Temp_class)23 0.004014 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lg_mss a.(T_) SexM   Densty
## log_mass    -0.574                            
## as.fc(T_)23 -0.512  0.532                     
## SexM        -0.196  0.216  0.080              
## Density     -0.768  0.109  0.107 -0.021       
## lg_:.(T_)23  0.402 -0.651 -0.759 -0.141 -0.111
ranova(MMR_mod_trtint)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log10(MMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (1 | Temp_class) + log_mass:as.factor(Temp_class)
##                  npar logLik     AIC        LRT Df Pr(>Chisq)
## <none>              8  534.2 -1052.4                         
## (1 | Temp_class)    7  534.2 -1054.4 6.8212e-13  1          1
# model with random slopes by temp treatment
MMR_mod_trt <- lmer(log10(MMR) ~ log_mass*as.factor(Temp_class) + Sex + Density + (0 + Temp_class|ID_fish),  
                 data=ds1, REML=T)

plot(MMR_mod_trt)

summary(MMR_mod_trt) 
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: log10(MMR) ~ log_mass * as.factor(Temp_class) + Sex + Density +  
##     (0 + Temp_class | ID_fish)
##    Data: ds1
## 
## REML criterion at convergence: -1086.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0147 -0.6251  0.0075  0.5777  3.5188 
## 
## Random effects:
##  Groups   Name         Variance  Std.Dev. Corr
##  ID_fish  Temp_class18 0.0008129 0.02851      
##           Temp_class23 0.0009152 0.03025  0.43
##  Residual              0.0042424 0.06513      
## Number of obs: 453, groups:  ID_fish, 91
## 
## Fixed effects:
##                                    Estimate Std. Error         df t value
## (Intercept)                       -0.034292   0.025461 208.686591  -1.347
## log_mass                           0.718885   0.022540 289.137042  31.894
## as.factor(Temp_class)23            0.051319   0.019279 235.372611   2.662
## SexM                               0.026620   0.009019  86.512834   2.952
## Density                           -0.012647   0.002880 209.059535  -4.391
## log_mass:as.factor(Temp_class)23  -0.097241   0.035385 315.337174  -2.748
##                                  Pr(>|t|)    
## (Intercept)                       0.17950    
## log_mass                          < 2e-16 ***
## as.factor(Temp_class)23           0.00831 ** 
## SexM                              0.00407 ** 
## Density                          1.79e-05 ***
## log_mass:as.factor(Temp_class)23  0.00634 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) lg_mss a.(T_) SexM   Densty
## log_mass    -0.604                            
## as.fc(T_)23 -0.451  0.609                     
## SexM        -0.204  0.180  0.062              
## Density     -0.839  0.151  0.128 -0.012       
## lg_:.(T_)23  0.396 -0.639 -0.883 -0.117 -0.109
ranova(MMR_mod_trt)
## ANOVA-like table for random-effects: Single term deletions
## 
## Model:
## log10(MMR) ~ log_mass + as.factor(Temp_class) + Sex + Density + (0 + Temp_class | ID_fish) + log_mass:as.factor(Temp_class)
##                                          npar logLik     AIC      LRT Df
## <none>                                     10 543.33 -1066.7            
## Temp_class in (0 + Temp_class | ID_fish)    8 543.31 -1070.6 0.035729  2
##                                          Pr(>Chisq)
## <none>                                             
## Temp_class in (0 + Temp_class | ID_fish)     0.9823
# test sex terms for model inclusion
MMR_mod_sex <- lmer(log10(MMR) ~ 0 + as.factor(Sex):as.factor(Temp_class) + as.factor(Sex):as.factor(Temp_class):log_mass + Density + (1 |ID_fish),  
                 data=ds1, REML=T)
plot(MMR_mod_sex)

summary(MMR_mod_sex)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## log10(MMR) ~ 0 + as.factor(Sex):as.factor(Temp_class) + as.factor(Sex):as.factor(Temp_class):log_mass +  
##     Density + (1 | ID_fish)
##    Data: ds1
## 
## REML criterion at convergence: -1073.9
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1037 -0.6119  0.0047  0.5638  3.5177 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev.
##  ID_fish  (Intercept) 0.0008985 0.02997 
##  Residual             0.0042347 0.06507 
## Number of obs: 453, groups:  ID_fish, 91
## 
## Fixed effects:
##                                                    Estimate Std. Error
## Density                                           -0.012730   0.002902
## as.factor(Sex)F:as.factor(Temp_class)18           -0.041773   0.027406
## as.factor(Sex)M:as.factor(Temp_class)18            0.008437   0.030355
## as.factor(Sex)F:as.factor(Temp_class)23            0.006331   0.029249
## as.factor(Sex)M:as.factor(Temp_class)23            0.048888   0.026927
## as.factor(Sex)F:as.factor(Temp_class)18:log_mass   0.731205   0.026963
## as.factor(Sex)M:as.factor(Temp_class)18:log_mass   0.688049   0.043257
## as.factor(Sex)F:as.factor(Temp_class)23:log_mass   0.650875   0.046066
## as.factor(Sex)M:as.factor(Temp_class)23:log_mass   0.606438   0.033600
##                                                          df t value Pr(>|t|)
## Density                                          211.516807  -4.386 1.82e-05
## as.factor(Sex)F:as.factor(Temp_class)18          209.572960  -1.524   0.1290
## as.factor(Sex)M:as.factor(Temp_class)18          281.936327   0.278   0.7812
## as.factor(Sex)F:as.factor(Temp_class)23          199.134068   0.216   0.8289
## as.factor(Sex)M:as.factor(Temp_class)23          206.145308   1.816   0.0709
## as.factor(Sex)F:as.factor(Temp_class)18:log_mass 357.870749  27.119  < 2e-16
## as.factor(Sex)M:as.factor(Temp_class)18:log_mass 439.391714  15.906  < 2e-16
## as.factor(Sex)F:as.factor(Temp_class)23:log_mass 292.475337  14.129  < 2e-16
## as.factor(Sex)M:as.factor(Temp_class)23:log_mass 365.600428  18.049  < 2e-16
##                                                     
## Density                                          ***
## as.factor(Sex)F:as.factor(Temp_class)18             
## as.factor(Sex)M:as.factor(Temp_class)18             
## as.factor(Sex)F:as.factor(Temp_class)23             
## as.factor(Sex)M:as.factor(Temp_class)23          .  
## as.factor(Sex)F:as.factor(Temp_class)18:log_mass ***
## as.factor(Sex)M:as.factor(Temp_class)18:log_mass ***
## as.factor(Sex)F:as.factor(Temp_class)23:log_mass ***
## as.factor(Sex)M:as.factor(Temp_class)23:log_mass ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##                 Densty as.(S)F:.(T_)18 as.(S)M:.(T_)18 as.(S)F:.(T_)23
## as.(S)F:.(T_)18 -0.790                                                
## as.(S)M:.(T_)18 -0.711  0.562                                         
## as.(S)F:.(T_)23 -0.620  0.490           0.440                         
## as.(S)M:.(T_)23 -0.731  0.578           0.520           0.453         
## a.(S)F:.(T_)18:  0.133 -0.656          -0.095          -0.083         
## a.(S)M:.(T_)18:  0.080 -0.063          -0.688          -0.050         
## a.(S)F:.(T_)23: -0.049  0.039           0.035          -0.682         
## a.(S)M:.(T_)23:  0.015 -0.012          -0.011          -0.009         
##                 as.(S)M:.(T_)23 a.(S)F:.(T_)18: a.(S)M:.(T_)18: a.(S)F:.(T_)23:
## as.(S)F:.(T_)18                                                                
## as.(S)M:.(T_)18                                                                
## as.(S)F:.(T_)23                                                                
## as.(S)M:.(T_)23                                                                
## a.(S)F:.(T_)18: -0.097                                                         
## a.(S)M:.(T_)18: -0.059           0.011                                         
## a.(S)F:.(T_)23:  0.036          -0.007          -0.004                         
## a.(S)M:.(T_)23: -0.580           0.002           0.001          -0.001
# plot MMR vs. log mass mean predicted level

pred_MMR <- ggpredict(MMR_mod, terms = c("log_mass", "Temp_class"))

ggplot() +
  geom_errorbar(data = pred_MMR, aes(x=x, ymin = conf.low, ymax = conf.high, color = group, group = group), width = 0.2) +   # add error bars showing 95% confidence intervals
  geom_point(data = pred_MMR, aes(x = x, y = predicted, color = group, group = group)) +  
  geom_line(data = pred_MMR, aes(x = x, y = predicted, color = group, group = group)) +
  labs(x = "log_mass", y = "Predicted MMR", color = "Temp_class") +
  theme_bw() +
  scale_color_manual(values = c("#78B7C5", "#F21A00")) +
  scale_y_log10() 

# Plot the raw data as a cloud around the predictions
ggplot() +
  geom_point(data = ds1, aes(x = log_mass, y = MMR, group = factor(Temp_class), color = factor(Temp_class)), 
             alpha = 0.3, position = position_jitter(width = 0.1)) +
  geom_point(data = pred_MMR, aes(x = x, y = predicted, color = group, group = group), 
             size = 3) +
  geom_line(data = pred_MMR, aes(x = x, y = predicted, color = group, group = group), 
            size = 2) +
  labs(x = "log_mass_centre", y = "Predicted MMR", color = "Temp_class") +
  theme_bw() +
  scale_color_manual(values = c("#78B7C5", "#F21A00")) +
  scale_y_log10() 

# F test on interaction effect

anova(MMR_mod, type = "3")
## Type III Analysis of Variance Table with Satterthwaite's method
##                                Sum Sq Mean Sq NumDF  DenDF   F value    Pr(>F)
## log_mass                       6.1674  6.1674     1 349.13 1453.5364 < 2.2e-16
## as.factor(Temp_class)          0.0296  0.0296     1 264.40    6.9784  0.008741
## Sex                            0.0367  0.0367     1  87.00    8.6421  0.004206
## Density                        0.0811  0.0811     1 211.71   19.1076 1.936e-05
## log_mass:as.factor(Temp_class) 0.0314  0.0314     1 359.55    7.4113  0.006797
##                                   
## log_mass                       ***
## as.factor(Temp_class)          ** 
## Sex                            ** 
## Density                        ***
## log_mass:as.factor(Temp_class) ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Generate fitted data

z_MMR<-augment(MMR_mod)                             # get variables and fitted values from model
head(z_MMR)
## # A tibble: 6 × 18
##   .rownames `log10(MMR)` log_mass `as.factor(Temp_class)` Sex   Density ID_fish
##   <chr>            <dbl>    <dbl> <fct>                   <fct>   <dbl>   <dbl>
## 1 1               0.154    0.423  23                      F           8      81
## 2 2               0.0458  -0.0915 23                      M           8      82
## 3 3               0.159    0.250  23                      F           8      83
## 4 4               0.215    0.436  23                      M           8      84
## 5 5               0.227    0.360  23                      M           8      85
## 6 7               0.0919   0.310  23                      M           8      87
## # ℹ 11 more variables: .fitted <dbl>, .resid <dbl>, .hat <dbl>, .cooksd <dbl>,
## #   .fixed <dbl>, .mu <dbl>, .offset <dbl>, .sqrtXwt <dbl>, .sqrtrwt <dbl>,
## #   .weights <dbl>, .wtres <dbl>
z_MMR <- z_MMR %>% 
  mutate(
    pwr.fitted = 10^(.fitted)
  )

# check assumptions: plot all residual values against the fitted 
ggplot(z_MMR, aes(x = .fitted, y = .resid)) +       # same plot as before
  geom_point(size = 2) 

ggplot(z_MMR, aes(x = log_mass, y = MMR, group = ID_fish)) + 
  geom_line(aes(y = pwr.fitted, color = `as.factor(Temp_class)`),            
            alpha = 0.6) +
  theme_bw() +
  scale_y_log10()+
  scale_color_manual(values  = c("#78B7C5", "#F21A00"))

FIGURE 2

Generate a figure that displays mass-scaling relationships for all three metabolic rates at both temperatures

### FIG. 2 ###

fig2.1 <- ggplot() +
    geom_ribbon(data = pred_SMR, aes(x=x, ymin = conf.low, ymax = conf.high, fill = group, color = group, group = group), width = 0.2, alpha = 0.1) +   # add error bars showing 95% confidence intervals
  geom_line(data = pred_SMR, aes(x = x, y = predicted, color = group, group = group), 
            size = 1) +
  geom_point(data = ds1, aes(x = log_mass, y = SMR, group = factor(Temp_class), color = factor(Temp_class)), 
             alpha = 0.25, position = position_jitter(width = 0.1)) +
  labs(x = "log10(mass (g))", y = expression(SMR~(mg~O[2]~h^-1)), color = expression("Temperature ("*degree*C*")")) +
  theme_classic() +
  scale_color_manual(values = c("#78B7C5", "#F21A00")) +
  scale_fill_manual(values = c("#78B7C5", "#F21A00"), guide = "none") +
  xlim(-0.25, 1.1) +
 scale_y_log10(limits = c(0.15, 1.6))

fig2.2 <- ggplot() +
    geom_ribbon(data = pred_RMR, aes(x=x, ymin = conf.low, ymax = conf.high, fill = group, color = group, group = group), width = 0.2, alpha = 0.1) +   # add error bars showing 95% confidence intervals
  geom_line(data = pred_RMR, aes(x = x, y = predicted, color = group, group = group), 
            size = 1) +
  geom_point(data = ds1, aes(x = log_mass, y = RMR, group = factor(Temp_class), color = factor(Temp_class)), 
             alpha = 0.25, position = position_jitter(width = 0.1)) +
  labs(x = "log10(mass (g))", y = expression(RMR~(mg~O[2]~h^-1)), color = expression("Temperature ("*degree*C*")")) +
  theme_classic() +
  scale_color_manual(values = c("#78B7C5", "#F21A00")) +
  scale_fill_manual(values = c("#78B7C5", "#F21A00"), guide = "none") +
  xlim(-0.25, 1.1) +
  scale_y_log10(limits=c(0.25,2.1))

fig2.3 <- ggplot() +
    geom_ribbon(data = pred_MMR, aes(x=x, ymin = conf.low, ymax = conf.high, fill = group, color = group, group = group), width = 0.2, alpha = 0.1) +   # add error bars showing 95% confidence intervals
  geom_line(data = pred_MMR, aes(x = x, y = predicted, color = group, group = group), 
            size = 1) +
  geom_point(data = ds1, aes(x = log_mass, y = MMR, group = factor(Temp_class), color = factor(Temp_class)), 
             alpha = 0.25, position = position_jitter(width = 0.1)) +
  labs(x = expression(log[10]~(mass~(g))), y = expression(MMR~(mg~O[2]~h^-1)), color = expression("Temperature ("*degree*C*")")) +
  theme_classic() +
  scale_color_manual(values = c("#78B7C5", "#F21A00")) +
  scale_fill_manual(values = c("#78B7C5", "#F21A00"), guide = "none") +
  xlim(-0.25, 1.1) +
  scale_y_log10(limits=c(0.5, 5))

fig2.4 <- ggplot(z_SMR, aes(x = log_mass, y = SMR, group = ID_fish)) + 
  geom_line(aes(y = pwr.fitted, color = `as.factor(Temp_class)`),      # note this will look odd if there are additional fixed effects      
            alpha = 0.6) +
  theme_classic() +
  scale_color_manual(values  = c("#78B7C5", "#F21A00")) +
  scale_y_log10(limits = c(0.15, 1.6))+
  xlim(-0.25, 1.1) +
  labs(x = expression(log[10]~(mass~(g))), y = expression(Fitted~SMR~(mg~O[2]~h^-1)), color = expression("Temperature ("*degree*C*")"))


fig2.5 <- ggplot(z_RMR, aes(x = log_mass, y = RMR, group = ID_fish)) + 
  geom_line(aes(y = pwr.fitted, color = `as.factor(Temp_class)`),      # note this will look odd if there are additional fixed effects      
            alpha = 0.6) +
  theme_classic() +
  scale_color_manual(values  = c("#78B7C5", "#F21A00")) +
  scale_y_log10(limits=c(0.25,2.1))+
  xlim(-0.25, 1.1) +
  labs(x = expression(log[10]~(mass~(g))), y = expression(Fitted~RMR~(mg~O[2]~h^-1)), color = expression("Temperature ("*degree*C*")"))

fig2.6 <- ggplot(z_MMR, aes(x = log_mass, y = MMR, group = ID_fish)) + 
  geom_line(aes(y = pwr.fitted, color = `as.factor(Temp_class)`),      # note this will look odd if there are additional fixed effects      
            alpha = 0.6) +
  theme_classic() +
  scale_color_manual(values  = c("#78B7C5", "#F21A00")) +
  scale_y_log10(limits=c(0.5, 5))+
  xlim(-0.25, 1.1) +
  labs(x = expression(log[10]~(mass~(g))), y = expression(Fitted~MMR~(mg~O[2]~h^-1)), color = expression("Temperature ("*degree*C*")"))

library(ggpubr)

fig2 <- ggarrange(print(fig2.1 + rremove("xlab")), 
                  print(fig2.2 + rremove("xlab")), 
                  print(fig2.3 + rremove("xlab")),
                  print(fig2.4 + rremove("xlab") + rremove("ylab")),
                  print(fig2.5 + rremove("xlab") + rremove("ylab")),
                  print(fig2.6 + rremove("xlab") + rremove("ylab")),
                           ncol = 3, nrow = 2, widths = c(5, 5, 5),
                           align = "v", 
                           common.legend = TRUE, 
                           legend = "right")

require(grid)
annotate_figure(fig2, bottom = textGrob(expression(log[10]~(mass~(g))), gp = gpar(cex = 1.1)))

setwd(figures_wd)
jpeg(file = "Fig2.jpeg", width = 1440, height=960, units = "px")
annotate_figure(fig2, bottom = textGrob(expression(log[10]~(mass~(g))), gp = gpar(cex = 1.1)))
dev.off()
## png 
##   2

UNIVARIATE MODEL OF FAT MASS

#subset ds1 to only include fish for which fat mass was measured
ds_fat <- ds1[which(!is.na(ds1$log_fatmass) & ds1$Month == "September"),]

fat_mod <- lm(log_fatmass ~ Sex + Sex:log_finalmass_centre + Sex:Temp_class + Temp_class + log_finalmass_centre:Temp_class + Sex:Temp_class:log_finalmass_centre + Density, data=ds_fat)

fat_mod1 <- lm(log_fatmass ~ Temp_class + log_finalmass_centre:Temp_class + SMR:Temp_class + Sex + Density, data=ds_fat) #using only final measurement of SMR here

fat_mod2 <- lm(log_fatmass ~ Temp_class + log_finalmass_centre:Temp_class + RMR:Temp_class + Sex + Density, data=ds_fat) #using only final measurement of RMR here

fat_mod3 <- lm(log_fatmass ~ Temp_class + log_finalmass_centre:Temp_class + MMR:Temp_class + Sex + Density, data=ds_fat) #using only final measurement of MMR here

summary(fat_mod)   # generate model output
## 
## Call:
## lm(formula = log_fatmass ~ Sex + Sex:log_finalmass_centre + Sex:Temp_class + 
##     Temp_class + log_finalmass_centre:Temp_class + Sex:Temp_class:log_finalmass_centre + 
##     Density, data = ds_fat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.66014 -0.10148  0.03335  0.10829  0.39864 
## 
## Coefficients:
##                                         Estimate Std. Error t value Pr(>|t|)
## (Intercept)                            -1.633537   0.145266 -11.245  < 2e-16
## SexM                                   -0.284274   0.179915  -1.580    0.118
## Temp_class23                           -0.207258   0.198303  -1.045    0.299
## Density                                -0.008147   0.013416  -0.607    0.545
## SexF:log_finalmass_centre               1.480710   0.186207   7.952 9.34e-12
## SexM:log_finalmass_centre               2.108583   0.246715   8.547 6.26e-13
## SexM:Temp_class23                      -0.124996   0.272443  -0.459    0.648
## log_finalmass_centre:Temp_class23       0.616740   0.357173   1.727    0.088
## SexM:log_finalmass_centre:Temp_class23 -0.030108   0.491759  -0.061    0.951
##                                           
## (Intercept)                            ***
## SexM                                      
## Temp_class23                              
## Density                                   
## SexF:log_finalmass_centre              ***
## SexM:log_finalmass_centre              ***
## SexM:Temp_class23                         
## log_finalmass_centre:Temp_class23      .  
## SexM:log_finalmass_centre:Temp_class23    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1789 on 81 degrees of freedom
##   (23 observations deleted due to missingness)
## Multiple R-squared:  0.807,  Adjusted R-squared:  0.7879 
## F-statistic: 42.34 on 8 and 81 DF,  p-value: < 2.2e-16
plot(fat_mod)      # residuals plot

anova(fat_mod)
## Analysis of Variance Table
## 
## Response: log_fatmass
##                                     Df  Sum Sq Mean Sq  F value  Pr(>F)    
## Sex                                  1  0.2138  0.2138   6.6797 0.01155 *  
## Temp_class                           1  0.0241  0.0241   0.7533 0.38801    
## Density                              1  0.0041  0.0041   0.1276 0.72191    
## Sex:log_finalmass_centre             2 10.3111  5.1556 161.0605 < 2e-16 ***
## Sex:Temp_class                       1  0.0995  0.0995   3.1072 0.08172 .  
## log_finalmass_centre:Temp_class      1  0.1894  0.1894   5.9183 0.01719 *  
## Sex:log_finalmass_centre:Temp_class  1  0.0001  0.0001   0.0037 0.95133    
## Residuals                           81  2.5928  0.0320                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#Marginal Effects
#Plotting marginal effects using the code Pete gave me
require(ggeffects)
pred.fat_mod <- ggpredict(fat_mod, terms = c("log_finalmass_centre", "Sex", "Temp_class"))
plot(pred.fat_mod)

ggplot() +
  geom_line(data = pred.fat_mod, aes(x = x, y = predicted, color = group, group = group)) +
  geom_ribbon(data = pred.fat_mod, aes(x= x, ymin = conf.low, ymax = conf.high, fill = group, group = group), alpha = 0.4) +
  labs(x = expression(log[10]~(centred~mass~(g))), y = expression(log[10]~(fat~mass~(g))), color = expression("Sex")) +
  theme_classic() +
  scale_color_manual(values = c("darkgreen", "purple")) +
  scale_fill_manual(values = c("darkgreen", "purple")) 

summary(fat_mod1)   # generate model output
## 
## Call:
## lm(formula = log_fatmass ~ Temp_class + log_finalmass_centre:Temp_class + 
##     SMR:Temp_class + Sex + Density, data = ds_fat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.75613 -0.09036  0.04041  0.12172  0.31448 
## 
## Coefficients:
##                                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                       -1.669290   0.132302 -12.617  < 2e-16 ***
## Temp_class23                      -0.335509   0.143164  -2.344   0.0215 *  
## SexM                              -0.028106   0.045131  -0.623   0.5352    
## Density                           -0.009403   0.014587  -0.645   0.5210    
## Temp_class18:log_finalmass_centre  1.373244   0.306634   4.478 2.41e-05 ***
## Temp_class23:log_finalmass_centre  2.827861   0.397690   7.111 3.91e-10 ***
## Temp_class18:SMR                   0.136718   0.146222   0.935   0.3525    
## Temp_class23:SMR                  -0.268150   0.236012  -1.136   0.2592    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1866 on 82 degrees of freedom
##   (23 observations deleted due to missingness)
## Multiple R-squared:  0.7876, Adjusted R-squared:  0.7694 
## F-statistic: 43.43 on 7 and 82 DF,  p-value: < 2.2e-16
plot(fat_mod)      # residuals plot

summary(fat_mod2)   # generate model output
## 
## Call:
## lm(formula = log_fatmass ~ Temp_class + log_finalmass_centre:Temp_class + 
##     RMR:Temp_class + Sex + Density, data = ds_fat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.76919 -0.08987  0.04334  0.12239  0.30783 
## 
## Coefficients:
##                                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                       -1.67164    0.13345 -12.526  < 2e-16 ***
## Temp_class23                      -0.33243    0.14187  -2.343   0.0215 *  
## SexM                              -0.02392    0.04425  -0.540   0.5904    
## Density                           -0.01111    0.01425  -0.780   0.4379    
## Temp_class18:log_finalmass_centre  1.40379    0.31899   4.401 3.22e-05 ***
## Temp_class23:log_finalmass_centre  2.74098    0.41355   6.628 3.35e-09 ***
## Temp_class18:RMR                   0.09706    0.12531   0.775   0.4409    
## Temp_class23:RMR                  -0.15521    0.18365  -0.845   0.4005    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1876 on 82 degrees of freedom
##   (23 observations deleted due to missingness)
## Multiple R-squared:  0.7853, Adjusted R-squared:  0.7669 
## F-statistic: 42.84 on 7 and 82 DF,  p-value: < 2.2e-16
plot(fat_mod)      # residuals plot

summary(fat_mod3)   # generate model output
## 
## Call:
## lm(formula = log_fatmass ~ Temp_class + log_finalmass_centre:Temp_class + 
##     MMR:Temp_class + Sex + Density, data = ds_fat)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.75804 -0.09095  0.02879  0.12219  0.33877 
## 
## Coefficients:
##                                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                       -1.666135   0.133622 -12.469  < 2e-16 ***
## Temp_class23                      -0.320143   0.141473  -2.263   0.0263 *  
## SexM                              -0.009131   0.042876  -0.213   0.8319    
## Density                           -0.014221   0.014104  -1.008   0.3163    
## Temp_class18:log_finalmass_centre  1.787541   0.296595   6.027 4.57e-08 ***
## Temp_class23:log_finalmass_centre  2.718573   0.301131   9.028 6.33e-14 ***
## Temp_class18:MMR                  -0.035603   0.058075  -0.613   0.5415    
## Temp_class23:MMR                  -0.075744   0.061506  -1.231   0.2217    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1869 on 82 degrees of freedom
##   (23 observations deleted due to missingness)
## Multiple R-squared:  0.7867, Adjusted R-squared:  0.7685 
## F-statistic: 43.21 on 7 and 82 DF,  p-value: < 2.2e-16
plot(fat_mod)      # residuals plot

pred_fat <- ggpredict(fat_mod, terms = c("log_finalmass_centre", "Temp_class"))

# Plot the raw data
ggplot() +
  geom_point(data = ds_fat, aes(x = log_finalmass_centre, y = log_fatmass, group = factor(Temp_class), color = factor(Temp_class))) +
  geom_line(data = pred_fat, aes(x = x, y = predicted, color = group, group = group)) +
  geom_ribbon(data = pred_fat, aes(x= x, ymin = conf.low, ymax = conf.high, fill = group, group = group), alpha = 0.4) +
  labs(x = expression(log[10]~(centred~mass~(g))), y = expression(log[10]~(fat~mass~(g))), color = expression("Temperature ("*degree*C*")")) +
  theme_classic() +
  scale_color_manual(values = c("#78B7C5", "#F21A00")) +
  scale_fill_manual(values = c("#78B7C5", "#F21A00")) 

ggplot(data = ds_fat, aes(y = Fat_mass_g, x = Sex, fill = factor(Sex))) +
  scale_fill_manual(values = c("darkgreen", "purple")) +
  geom_boxplot() +
  labs(x = "Sex", y = "Fat mass (g)") +
  theme_classic()

# F test on interaction term
anova(fat_mod)
## Analysis of Variance Table
## 
## Response: log_fatmass
##                                     Df  Sum Sq Mean Sq  F value  Pr(>F)    
## Sex                                  1  0.2138  0.2138   6.6797 0.01155 *  
## Temp_class                           1  0.0241  0.0241   0.7533 0.38801    
## Density                              1  0.0041  0.0041   0.1276 0.72191    
## Sex:log_finalmass_centre             2 10.3111  5.1556 161.0605 < 2e-16 ***
## Sex:Temp_class                       1  0.0995  0.0995   3.1072 0.08172 .  
## log_finalmass_centre:Temp_class      1  0.1894  0.1894   5.9183 0.01719 *  
## Sex:log_finalmass_centre:Temp_class  1  0.0001  0.0001   0.0037 0.95133    
## Residuals                           81  2.5928  0.0320                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(fat_mod1)
## Analysis of Variance Table
## 
## Response: log_fatmass
##                                 Df  Sum Sq Mean Sq  F value  Pr(>F)    
## Temp_class                       1  0.0406  0.0406   1.1650 0.28359    
## Sex                              1  0.1974  0.1974   5.6706 0.01957 *  
## Density                          1  0.0041  0.0041   0.1173 0.73285    
## Temp_class:log_finalmass_centre  2 10.2619  5.1309 147.4086 < 2e-16 ***
## Temp_class:SMR                   2  0.0769  0.0384   1.1045 0.33625    
## Residuals                       82  2.8542  0.0348                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(fat_mod2)
## Analysis of Variance Table
## 
## Response: log_fatmass
##                                 Df  Sum Sq Mean Sq  F value  Pr(>F)    
## Temp_class                       1  0.0406  0.0406   1.1526 0.28616    
## Sex                              1  0.1974  0.1974   5.6100 0.02021 *  
## Density                          1  0.0041  0.0041   0.1161 0.73423    
## Temp_class:log_finalmass_centre  2 10.2619  5.1309 145.8320 < 2e-16 ***
## Temp_class:RMR                   2  0.0460  0.0230   0.6542 0.52257    
## Residuals                       82  2.8851  0.0352                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
anova(fat_mod3)
## Analysis of Variance Table
## 
## Response: log_fatmass
##                                 Df  Sum Sq Mean Sq  F value  Pr(>F)    
## Temp_class                       1  0.0406  0.0406   1.1604 0.28455    
## Sex                              1  0.1974  0.1974   5.6480 0.01981 *  
## Density                          1  0.0041  0.0041   0.1168 0.73337    
## Temp_class:log_finalmass_centre  2 10.2619  5.1309 146.8201 < 2e-16 ***
## Temp_class:MMR                   2  0.0654  0.0327   0.9364 0.39619    
## Residuals                       82  2.8657  0.0349                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

This model is excluded from the following covariance matrix because we were only able to collect one measurement of fat mass from each fish, all at the very end of the experiment. Without an initial measurement of fat mass (logistically impossible as it was a lethal measurement), we are unable to assess within-subject residual correlation. For this reason, we have assessed it separately.

COVARIANCE MATRIX

covariance matrix with only sex specific traits

RN <- bf(mass1 ~ 0 + Sex:Temp_class + Sex:Temp_class:Month_centred_factor + Density + (1 + Month_centred |a| ID_fish)) + 
      bf(SMR1  ~ 0 + Sex:Temp_class + Sex:Temp_class:log_mass_centre + Density + (1 |a| ID_fish)) + 
      bf(RMR1  ~ 0 + Sex:Temp_class + Sex:Temp_class:log_mass_centre + Density + (1 |a| ID_fish)) + 
      bf(MMR1  ~ 0 + Sex:Temp_class + Sex:Temp_class:log_mass_centre + Density + (1 |a| ID_fish)) + 
      set_rescor(T)

prior1 <- c(set_prior("normal(0,2)", class = "b", resp = c("mass1", "SMR1", "RMR1","MMR1")), 
            set_prior("cauchy(0,2)", class = "sd", resp = c("mass1", "SMR1", "RMR1","MMR1")), 
            set_prior("lkj(2)", class = "cor"))

 fit0 <- brm(formula = RN, 
            data = ds1, 
            prior = prior1, 
            seed = 42, warmup = 500, iter = 6000, chains = 4, cores = 4,
            control=list(adapt_delta=0.97, max_treedepth = 15),
            save_ranef = T
           )
 
print(fit0, digits = 3)
##  Family: MV(gaussian, gaussian, gaussian, gaussian) 
##   Links: mu = identity; sigma = identity
##          mu = identity; sigma = identity
##          mu = identity; sigma = identity
##          mu = identity; sigma = identity 
## Formula: mass1 ~ 0 + Sex:Temp_class + Sex:Temp_class:Month_centred_factor + Density + (1 + Month_centred | a | ID_fish) 
##          SMR1 ~ 0 + Sex:Temp_class + Sex:Temp_class:log_mass_centre + Density + (1 | a | ID_fish) 
##          RMR1 ~ 0 + Sex:Temp_class + Sex:Temp_class:log_mass_centre + Density + (1 | a | ID_fish) 
##          MMR1 ~ 0 + Sex:Temp_class + Sex:Temp_class:log_mass_centre + Density + (1 | a | ID_fish) 
##    Data: ds1 (Number of observations: 453) 
##   Draws: 4 chains, each with iter = 6000; warmup = 500; thin = 1;
##          total post-warmup draws = 22000
## 
## Multilevel Hyperparameters:
## ~ID_fish (Number of levels: 91) 
##                                          Estimate Est.Error l-95% CI u-95% CI
## sd(mass1_Intercept)                         0.167     0.014    0.142    0.197
## sd(mass1_Month_centred)                     0.039     0.004    0.032    0.047
## sd(SMR1_Intercept)                          0.012     0.007    0.001    0.026
## sd(RMR1_Intercept)                          0.032     0.005    0.023    0.043
## sd(MMR1_Intercept)                          0.031     0.005    0.021    0.041
## cor(mass1_Intercept,mass1_Month_centred)   -0.358     0.103   -0.545   -0.145
## cor(mass1_Intercept,SMR1_Intercept)        -0.110     0.304   -0.667    0.512
## cor(mass1_Month_centred,SMR1_Intercept)     0.331     0.286   -0.341    0.787
## cor(mass1_Intercept,RMR1_Intercept)         0.160     0.152   -0.149    0.444
## cor(mass1_Month_centred,RMR1_Intercept)     0.179     0.151   -0.121    0.464
## cor(SMR1_Intercept,RMR1_Intercept)          0.394     0.316   -0.375    0.832
## cor(mass1_Intercept,MMR1_Intercept)         0.179     0.167   -0.161    0.488
## cor(mass1_Month_centred,MMR1_Intercept)     0.199     0.153   -0.111    0.488
## cor(SMR1_Intercept,MMR1_Intercept)          0.008     0.308   -0.595    0.598
## cor(RMR1_Intercept,MMR1_Intercept)          0.399     0.165    0.051    0.694
##                                           Rhat Bulk_ESS Tail_ESS
## sd(mass1_Intercept)                      1.000     5430     9939
## sd(mass1_Month_centred)                  1.000     5852     9991
## sd(SMR1_Intercept)                       1.001     3809     8610
## sd(RMR1_Intercept)                       1.001     2154     6458
## sd(MMR1_Intercept)                       1.000     7495    11259
## cor(mass1_Intercept,mass1_Month_centred) 1.000     6044    10623
## cor(mass1_Intercept,SMR1_Intercept)      1.000    15888    15169
## cor(mass1_Month_centred,SMR1_Intercept)  1.000    17212    14454
## cor(mass1_Intercept,RMR1_Intercept)      1.000    11506    14354
## cor(mass1_Month_centred,RMR1_Intercept)  1.001    10437    15802
## cor(SMR1_Intercept,RMR1_Intercept)       1.003     1116     2049
## cor(mass1_Intercept,MMR1_Intercept)      1.000    12958    15452
## cor(mass1_Month_centred,MMR1_Intercept)  1.000    13864    16602
## cor(SMR1_Intercept,MMR1_Intercept)       1.001     2043     4458
## cor(RMR1_Intercept,MMR1_Intercept)       1.000     7519    14396
## 
## Regression Coefficients:
##                                               Estimate Est.Error l-95% CI
## mass1_Density                                    0.012     0.006    0.000
## mass1_SexF:Temp_class18                          0.284     0.054    0.180
## mass1_SexM:Temp_class18                          0.206     0.060    0.088
## mass1_SexF:Temp_class23                          0.230     0.057    0.119
## mass1_SexM:Temp_class23                          0.148     0.060    0.031
## mass1_SexF:Temp_class18:Month_centred_factor1    0.161     0.016    0.130
## mass1_SexM:Temp_class18:Month_centred_factor1    0.157     0.021    0.117
## mass1_SexF:Temp_class23:Month_centred_factor1    0.104     0.021    0.062
## mass1_SexM:Temp_class23:Month_centred_factor1    0.161     0.023    0.116
## mass1_SexF:Temp_class18:Month_centred_factor2    0.254     0.020    0.214
## mass1_SexM:Temp_class18:Month_centred_factor2    0.174     0.027    0.121
## mass1_SexF:Temp_class23:Month_centred_factor2    0.164     0.026    0.114
## mass1_SexM:Temp_class23:Month_centred_factor2    0.270     0.028    0.215
## mass1_SexF:Temp_class18:Month_centred_factor3    0.198     0.026    0.147
## mass1_SexM:Temp_class18:Month_centred_factor3    0.114     0.034    0.047
## mass1_SexF:Temp_class23:Month_centred_factor3    0.169     0.033    0.105
## mass1_SexM:Temp_class23:Month_centred_factor3    0.277     0.036    0.207
## mass1_SexF:Temp_class18:Month_centred_factor4    0.407     0.031    0.346
## mass1_SexM:Temp_class18:Month_centred_factor4    0.340     0.040    0.260
## mass1_SexF:Temp_class23:Month_centred_factor4    0.288     0.040    0.209
## mass1_SexM:Temp_class23:Month_centred_factor4    0.425     0.044    0.339
## SMR1_Density                                    -0.002     0.003   -0.009
## SMR1_SexF:Temp_class18                          -0.875     0.040   -0.953
## SMR1_SexM:Temp_class18                          -0.967     0.048   -1.062
## SMR1_SexF:Temp_class23                          -0.683     0.045   -0.772
## SMR1_SexM:Temp_class23                          -0.647     0.040   -0.726
## SMR1_SexF:Temp_class18:log_mass_centre           0.806     0.037    0.734
## SMR1_SexM:Temp_class18:log_mass_centre           0.963     0.058    0.850
## SMR1_SexF:Temp_class23:log_mass_centre           0.710     0.058    0.599
## SMR1_SexM:Temp_class23:log_mass_centre           0.674     0.045    0.584
## RMR1_Density                                    -0.006     0.003   -0.012
## RMR1_SexF:Temp_class18                          -0.615     0.039   -0.691
## RMR1_SexM:Temp_class18                          -0.689     0.045   -0.779
## RMR1_SexF:Temp_class23                          -0.466     0.045   -0.554
## RMR1_SexM:Temp_class23                          -0.436     0.039   -0.513
## RMR1_SexF:Temp_class18:log_mass_centre           0.710     0.035    0.642
## RMR1_SexM:Temp_class18:log_mass_centre           0.821     0.053    0.718
## RMR1_SexF:Temp_class23:log_mass_centre           0.622     0.056    0.514
## RMR1_SexM:Temp_class23:log_mass_centre           0.579     0.042    0.494
## MMR1_Density                                    -0.013     0.003   -0.019
## MMR1_SexF:Temp_class18                          -0.204     0.035   -0.272
## MMR1_SexM:Temp_class18                          -0.148     0.040   -0.225
## MMR1_SexF:Temp_class23                          -0.132     0.040   -0.210
## MMR1_SexM:Temp_class23                          -0.085     0.035   -0.153
## MMR1_SexF:Temp_class18:log_mass_centre           0.718     0.031    0.657
## MMR1_SexM:Temp_class18:log_mass_centre           0.678     0.045    0.590
## MMR1_SexF:Temp_class23:log_mass_centre           0.629     0.049    0.533
## MMR1_SexM:Temp_class23:log_mass_centre           0.593     0.037    0.519
##                                               u-95% CI  Rhat Bulk_ESS Tail_ESS
## mass1_Density                                    0.024 1.000     6711    10981
## mass1_SexF:Temp_class18                          0.390 1.000     4564     8643
## mass1_SexM:Temp_class18                          0.324 1.000     4386     8415
## mass1_SexF:Temp_class23                          0.343 1.000     4265     7980
## mass1_SexM:Temp_class23                          0.268 1.000     4726     7882
## mass1_SexF:Temp_class18:Month_centred_factor1    0.193 1.001    11148    14768
## mass1_SexM:Temp_class18:Month_centred_factor1    0.198 1.000    10838    15271
## mass1_SexF:Temp_class23:Month_centred_factor1    0.145 1.000    11252    14898
## mass1_SexM:Temp_class23:Month_centred_factor1    0.205 1.000    10625    14688
## mass1_SexF:Temp_class18:Month_centred_factor2    0.294 1.001     5942    11371
## mass1_SexM:Temp_class18:Month_centred_factor2    0.227 1.002     5705    10164
## mass1_SexF:Temp_class23:Month_centred_factor2    0.216 1.000     6792    11809
## mass1_SexM:Temp_class23:Month_centred_factor2    0.326 1.000     6631     9866
## mass1_SexF:Temp_class18:Month_centred_factor3    0.248 1.001     4768     8353
## mass1_SexM:Temp_class18:Month_centred_factor3    0.181 1.002     4628     9215
## mass1_SexF:Temp_class23:Month_centred_factor3    0.235 1.000     5556     9519
## mass1_SexM:Temp_class23:Month_centred_factor3    0.348 1.000     5297     8763
## mass1_SexF:Temp_class18:Month_centred_factor4    0.468 1.001     4489     8492
## mass1_SexM:Temp_class18:Month_centred_factor4    0.419 1.002     4337     8685
## mass1_SexF:Temp_class23:Month_centred_factor4    0.367 1.000     5037     9083
## mass1_SexM:Temp_class23:Month_centred_factor4    0.510 1.000     4961     7621
## SMR1_Density                                     0.004 1.000    12077    16800
## SMR1_SexF:Temp_class18                          -0.797 1.001     7008    11864
## SMR1_SexM:Temp_class18                          -0.874 1.001     7995    12588
## SMR1_SexF:Temp_class23                          -0.595 1.001     8792    12493
## SMR1_SexM:Temp_class23                          -0.570 1.001     7953    13161
## SMR1_SexF:Temp_class18:log_mass_centre           0.878 1.001     8564    13224
## SMR1_SexM:Temp_class18:log_mass_centre           1.076 1.000     9974    14180
## SMR1_SexF:Temp_class23:log_mass_centre           0.824 1.001    10179    14311
## SMR1_SexM:Temp_class23:log_mass_centre           0.760 1.000     9639    14921
## RMR1_Density                                     0.001 1.000    10974    15126
## RMR1_SexF:Temp_class18                          -0.537 1.001     7075    12036
## RMR1_SexM:Temp_class18                          -0.601 1.001     7864    12329
## RMR1_SexF:Temp_class23                          -0.379 1.001     8311    11983
## RMR1_SexM:Temp_class23                          -0.359 1.000     7899    13195
## RMR1_SexF:Temp_class18:log_mass_centre           0.778 1.001     8538    13846
## RMR1_SexM:Temp_class18:log_mass_centre           0.927 1.000    10273    13978
## RMR1_SexF:Temp_class23:log_mass_centre           0.733 1.000     9957    13423
## RMR1_SexM:Temp_class23:log_mass_centre           0.661 1.000     9832    14947
## MMR1_Density                                    -0.007 1.000    14150    16707
## MMR1_SexF:Temp_class18                          -0.135 1.000     9836    14603
## MMR1_SexM:Temp_class18                          -0.070 1.000    11607    14923
## MMR1_SexF:Temp_class23                          -0.054 1.000    11757    14948
## MMR1_SexM:Temp_class23                          -0.016 1.000    11091    15576
## MMR1_SexF:Temp_class18:log_mass_centre           0.779 1.000    12482    15761
## MMR1_SexM:Temp_class18:log_mass_centre           0.765 1.000    16329    16198
## MMR1_SexF:Temp_class23:log_mass_centre           0.725 1.000    15410    15767
## MMR1_SexM:Temp_class23:log_mass_centre           0.665 1.000    16232    16397
## 
## Further Distributional Parameters:
##             Estimate Est.Error l-95% CI u-95% CI  Rhat Bulk_ESS Tail_ESS
## sigma_mass1    0.060     0.003    0.055    0.065 1.000    13626    15029
## sigma_SMR1     0.091     0.003    0.085    0.097 1.000    15574    16610
## sigma_RMR1     0.080     0.003    0.075    0.087 1.000    12035    17014
## sigma_MMR1     0.066     0.003    0.061    0.071 1.000    18093    16365
## 
## Residual Correlations: 
##                    Estimate Est.Error l-95% CI u-95% CI  Rhat Bulk_ESS Tail_ESS
## rescor(mass1,SMR1)   -0.086     0.079   -0.239    0.070 1.000     8393    13542
## rescor(mass1,RMR1)   -0.108     0.074   -0.252    0.037 1.000     8805    13667
## rescor(SMR1,RMR1)     0.862     0.013    0.834    0.886 1.000    16325    17291
## rescor(mass1,MMR1)   -0.050     0.069   -0.185    0.086 1.000    14591    16962
## rescor(SMR1,MMR1)     0.143     0.050    0.046    0.239 1.001    26459    16631
## rescor(RMR1,MMR1)     0.207     0.050    0.108    0.303 1.001    26753    17238
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(fit0)

conditional_effects(fit0)

build covariance matrix with brms

RN <- bf(mass1 ~ 0 + Temp_class + Month_centred_factor:Temp_class + Sex + Density + (1 + Month_centred |a| ID_fish)) + 
      bf(SMR1  ~ 0 + Temp_class + log_mass_centre:Temp_class + Sex + Density + (1 |a| ID_fish)) + 
      bf(RMR1  ~ 0 + Temp_class + log_mass_centre:Temp_class + Sex + Density + (1 |a| ID_fish)) + 
      bf(MMR1  ~ 0 + Temp_class + log_mass_centre:Temp_class + Sex + Density + (1 |a| ID_fish)) + 
      set_rescor(T)

prior1 <- c(set_prior("normal(0,2)", class = "b", resp = c("mass1", "SMR1", "RMR1","MMR1")), 
            set_prior("cauchy(0,2)", class = "sd", resp = c("mass1", "SMR1", "RMR1","MMR1")), 
            set_prior("lkj(2)", class = "cor"))

 fit1 <- brm(formula = RN, 
            data = ds1, 
            prior = prior1, 
            seed = 42, warmup = 500, iter = 10000, chains = 4, cores = 4,
            control=list(adapt_delta=0.97, max_treedepth = 15),
            save_ranef = T
           )

print(fit1, digits = 3)
##  Family: MV(gaussian, gaussian, gaussian, gaussian) 
##   Links: mu = identity; sigma = identity
##          mu = identity; sigma = identity
##          mu = identity; sigma = identity
##          mu = identity; sigma = identity 
## Formula: mass1 ~ 0 + Temp_class + Month_centred_factor:Temp_class + Sex + Density + (1 + Month_centred | a | ID_fish) 
##          SMR1 ~ 0 + Temp_class + log_mass_centre:Temp_class + Sex + Density + (1 | a | ID_fish) 
##          RMR1 ~ 0 + Temp_class + log_mass_centre:Temp_class + Sex + Density + (1 | a | ID_fish) 
##          MMR1 ~ 0 + Temp_class + log_mass_centre:Temp_class + Sex + Density + (1 | a | ID_fish) 
##    Data: ds1 (Number of observations: 453) 
##   Draws: 4 chains, each with iter = 10000; warmup = 500; thin = 1;
##          total post-warmup draws = 38000
## 
## Multilevel Hyperparameters:
## ~ID_fish (Number of levels: 91) 
##                                          Estimate Est.Error l-95% CI u-95% CI
## sd(mass1_Intercept)                         0.165     0.014    0.141    0.194
## sd(mass1_Month_centred)                     0.040     0.004    0.034    0.048
## sd(SMR1_Intercept)                          0.011     0.007    0.001    0.025
## sd(RMR1_Intercept)                          0.031     0.005    0.023    0.042
## sd(MMR1_Intercept)                          0.030     0.005    0.020    0.040
## cor(mass1_Intercept,mass1_Month_centred)   -0.318     0.106   -0.512   -0.101
## cor(mass1_Intercept,SMR1_Intercept)        -0.128     0.307   -0.688    0.503
## cor(mass1_Month_centred,SMR1_Intercept)     0.338     0.293   -0.346    0.800
## cor(mass1_Intercept,RMR1_Intercept)         0.158     0.151   -0.148    0.440
## cor(mass1_Month_centred,RMR1_Intercept)     0.177     0.149   -0.117    0.466
## cor(SMR1_Intercept,RMR1_Intercept)          0.363     0.330   -0.428    0.825
## cor(mass1_Intercept,MMR1_Intercept)         0.191     0.169   -0.156    0.505
## cor(mass1_Month_centred,MMR1_Intercept)     0.169     0.153   -0.137    0.460
## cor(SMR1_Intercept,MMR1_Intercept)          0.017     0.311   -0.592    0.610
## cor(RMR1_Intercept,MMR1_Intercept)          0.411     0.165    0.063    0.705
##                                           Rhat Bulk_ESS Tail_ESS
## sd(mass1_Intercept)                      1.001     8724    15843
## sd(mass1_Month_centred)                  1.000    11231    18796
## sd(SMR1_Intercept)                       1.001     7027    11325
## sd(RMR1_Intercept)                       1.001     3711    10259
## sd(MMR1_Intercept)                       1.000    13870    16558
## cor(mass1_Intercept,mass1_Month_centred) 1.001     9327    18688
## cor(mass1_Intercept,SMR1_Intercept)      1.000    33502    27685
## cor(mass1_Month_centred,SMR1_Intercept)  1.000    28395    23325
## cor(mass1_Intercept,RMR1_Intercept)      1.000    21696    25906
## cor(mass1_Month_centred,RMR1_Intercept)  1.000    16783    23670
## cor(SMR1_Intercept,RMR1_Intercept)       1.001     1755     3518
## cor(mass1_Intercept,MMR1_Intercept)      1.000    23420    26235
## cor(mass1_Month_centred,MMR1_Intercept)  1.000    26510    28843
## cor(SMR1_Intercept,MMR1_Intercept)       1.000     3062     7410
## cor(RMR1_Intercept,MMR1_Intercept)       1.000    13833    24495
## 
## Regression Coefficients:
##                                          Estimate Est.Error l-95% CI u-95% CI
## mass1_Temp_class18                          0.290     0.052    0.188    0.392
## mass1_Temp_class23                          0.232     0.052    0.130    0.336
## mass1_SexM                                 -0.072     0.034   -0.139   -0.004
## mass1_Density                               0.011     0.006   -0.001    0.023
## mass1_Temp_class18:Month_centred_factor1    0.158     0.013    0.133    0.184
## mass1_Temp_class23:Month_centred_factor1    0.129     0.016    0.098    0.161
## mass1_Temp_class18:Month_centred_factor2    0.228     0.017    0.195    0.260
## mass1_Temp_class23:Month_centred_factor2    0.216     0.020    0.177    0.255
## mass1_Temp_class18:Month_centred_factor3    0.170     0.021    0.128    0.212
## mass1_Temp_class23:Month_centred_factor3    0.224     0.025    0.174    0.273
## mass1_Temp_class18:Month_centred_factor4    0.380     0.026    0.330    0.430
## mass1_Temp_class23:Month_centred_factor4    0.354     0.031    0.293    0.415
## SMR1_Temp_class18                          -0.897     0.036   -0.969   -0.826
## SMR1_Temp_class23                          -0.664     0.034   -0.731   -0.597
## SMR1_SexM                                   0.015     0.009   -0.004    0.034
## SMR1_Density                               -0.002     0.003   -0.009    0.004
## SMR1_Temp_class18:log_mass_centre           0.836     0.032    0.772    0.899
## SMR1_Temp_class23:log_mass_centre           0.681     0.037    0.608    0.753
## RMR1_Temp_class18                          -0.631     0.036   -0.701   -0.561
## RMR1_Temp_class23                          -0.446     0.034   -0.512   -0.378
## RMR1_SexM                                   0.002     0.010   -0.018    0.023
## RMR1_Density                               -0.006     0.003   -0.012    0.001
## RMR1_Temp_class18:log_mass_centre           0.732     0.030    0.673    0.791
## RMR1_Temp_class23:log_mass_centre           0.592     0.035    0.523    0.662
## MMR1_Temp_class18                          -0.195     0.033   -0.259   -0.131
## MMR1_Temp_class23                          -0.121     0.031   -0.182   -0.060
## MMR1_SexM                                   0.026     0.009    0.008    0.044
## MMR1_Density                               -0.013     0.003   -0.019   -0.007
## MMR1_Temp_class18:log_mass_centre           0.708     0.027    0.654    0.762
## MMR1_Temp_class23:log_mass_centre           0.609     0.032    0.547    0.672
##                                           Rhat Bulk_ESS Tail_ESS
## mass1_Temp_class18                       1.000     9281    16767
## mass1_Temp_class23                       1.000     8639    15827
## mass1_SexM                               1.000     5903    10692
## mass1_Density                            1.000    14045    22047
## mass1_Temp_class18:Month_centred_factor1 1.000    18640    26859
## mass1_Temp_class23:Month_centred_factor1 1.000    18861    25964
## mass1_Temp_class18:Month_centred_factor2 1.000    10533    19970
## mass1_Temp_class23:Month_centred_factor2 1.000    11106    20122
## mass1_Temp_class18:Month_centred_factor3 1.000     8770    15803
## mass1_Temp_class23:Month_centred_factor3 1.000     9191    16982
## mass1_Temp_class18:Month_centred_factor4 1.000     8096    15500
## mass1_Temp_class23:Month_centred_factor4 1.000     8171    14915
## SMR1_Temp_class18                        1.000    16774    22576
## SMR1_Temp_class23                        1.000    16149    23956
## SMR1_SexM                                1.000    45285    29965
## SMR1_Density                             1.000    32156    28146
## SMR1_Temp_class18:log_mass_centre        1.000    19146    25494
## SMR1_Temp_class23:log_mass_centre        1.000    19297    25087
## RMR1_Temp_class18                        1.000    14900    21747
## RMR1_Temp_class23                        1.000    14417    22655
## RMR1_SexM                                1.000    26137    29528
## RMR1_Density                             1.000    25322    28582
## RMR1_Temp_class18:log_mass_centre        1.000    18250    25032
## RMR1_Temp_class23:log_mass_centre        1.000    17444    23199
## MMR1_Temp_class18                        1.000    19049    24584
## MMR1_Temp_class23                        1.000    17889    24390
## MMR1_SexM                                1.000    30258    30267
## MMR1_Density                             1.000    29008    30210
## MMR1_Temp_class18:log_mass_centre        1.000    23182    26563
## MMR1_Temp_class23:log_mass_centre        1.000    24019    27582
## 
## Further Distributional Parameters:
##             Estimate Est.Error l-95% CI u-95% CI  Rhat Bulk_ESS Tail_ESS
## sigma_mass1    0.060     0.003    0.055    0.066 1.000    22767    28921
## sigma_SMR1     0.091     0.003    0.085    0.098 1.000    24535    28184
## sigma_RMR1     0.081     0.003    0.075    0.087 1.000    18866    28313
## sigma_MMR1     0.066     0.002    0.061    0.071 1.000    30665    26391
## 
## Residual Correlations: 
##                    Estimate Est.Error l-95% CI u-95% CI  Rhat Bulk_ESS Tail_ESS
## rescor(mass1,SMR1)   -0.017     0.076   -0.166    0.133 1.000    17064    24877
## rescor(mass1,RMR1)   -0.056     0.073   -0.198    0.087 1.000    17052    25013
## rescor(SMR1,RMR1)     0.863     0.013    0.836    0.887 1.000    27400    30731
## rescor(mass1,MMR1)   -0.069     0.069   -0.203    0.065 1.000    26677    30302
## rescor(SMR1,MMR1)     0.138     0.049    0.040    0.234 1.000    45101    29153
## rescor(RMR1,MMR1)     0.203     0.049    0.105    0.298 1.000    45295    29787
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
plot(fit1)

conditional_effects(fit1)

REPEATABILITY ESTIMATES

Use posterior distribution to generate repeatability estimates

View(posterior::summarise_draws(as_draws_array(fit1))) # view all predicted values of random effects

# random effects are fitted as standard dev, not variance, so we need to square them first
# var names to insert below: SMR, RMR, MMR, growth

# SMR repeatability

Vint_SMR <- (as_draws_array(fit1, variable = "sd_ID_fish__SMR1_Intercept"))^2
Vresid_SMR <- (as_draws_array(fit1, variable = "sigma_SMR1"))^2
R_SMR <- Vint_SMR / (Vint_SMR + Vresid_SMR)

"SMR"; mean(R_SMR); quantile(R_SMR, probs = c(0.025, 0.975))
## [1] "SMR"
## [1] 0.01979299
##         2.5%        97.5% 
## 5.327465e-05 7.128585e-02
# RMR repeatability

Vint_RMR <- (as_draws_array(fit1, variable = "sd_ID_fish__RMR1_Intercept"))^2
Vresid_RMR <- (as_draws_array(fit1, variable = "sigma_RMR1"))^2
R_RMR <- Vint_RMR / (Vint_RMR + Vresid_RMR)

"RMR"; mean(R_RMR); quantile(R_RMR, probs = c(0.025, 0.975))
## [1] "RMR"
## [1] 0.130803
##       2.5%      97.5% 
## 0.06931255 0.21816931
# MMR repeatability

Vint_MMR <- (as_draws_array(fit1, variable = "sd_ID_fish__MMR1_Intercept"))^2
Vresid_MMR <- (as_draws_array(fit1, variable = "sigma_MMR1"))^2
R_MMR <- Vint_MMR / (Vint_MMR + Vresid_MMR)

"MMR"; mean(R_MMR); quantile(R_MMR, probs = c(0.025, 0.975))
## [1] "MMR"
## [1] 0.1705791
##       2.5%      97.5% 
## 0.07647187 0.27647865
# Growth rate repeatability
Vint_growth <- (as_draws_array(fit1, variable = "sd_ID_fish__mass1_Month_centred"))^2
Vslope_growth <- (as_draws_array(fit1, variable = "sd_ID_fish__mass1_Intercept"))^2
corr_growth <- (as_draws_array(fit1, variable = "cor_ID_fish__mass1_Intercept__mass1_Month_centred"))
Vresid_growth <- (as_draws_array(fit1, variable = "sigma_mass1"))^2
COVis_growth <- corr_growth*sqrt(Vslope_growth)*sqrt(Vint_growth)

x<- 0 #(can be 0 - 4)
VAR0 <- Vint_growth +2*COVis_growth*x + Vslope_growth*(x^2)
R0_growth <- VAR0/(VAR0 + Vresid_growth)

"Growth - May"; mean(R0_growth); quantile(R0_growth, probs = c(0.025, 0.975))
## [1] "Growth - May"
## [1] 0.3126487
##      2.5%     97.5% 
## 0.2276433 0.4086757
x<- 1 #(can be 0 - 4)
VAR1 <- Vint_growth +2*COVis_growth*x + Vslope_growth*(x^2)
R1_growth <- VAR1/(VAR1 + Vresid_growth)

"Growth - June"; mean(R1_growth); quantile(R1_growth, probs = c(0.025, 0.975))
## [1] "Growth - June"
## [1] 0.8700228
##      2.5%     97.5% 
## 0.8252644 0.9079288
x<- 2 #(can be 0 - 4)
VAR2 <- Vint_growth +2*COVis_growth*x + Vslope_growth*(x^2)
R2_growth <- VAR2/(VAR2 + Vresid_growth)

"Growth - July"; mean(R2_growth); quantile(R2_growth, probs = c(0.025, 0.975))
## [1] "Growth - July"
## [1] 0.9651094
##      2.5%     97.5% 
## 0.9511292 0.9762125
x<- 3 #(can be 0 - 4)
VAR3 <- Vint_growth +2*COVis_growth*x + Vslope_growth*(x^2)
R3_growth <- VAR3/(VAR3 + Vresid_growth)

"Growth - August"; mean(R3_growth); quantile(R3_growth, probs = c(0.025, 0.975))
## [1] "Growth - August"
## [1] 0.9844639
##      2.5%     97.5% 
## 0.9780419 0.9894882
x<- 4 #(can be 0 - 4)
VAR4 <- Vint_growth +2*COVis_growth*x + Vslope_growth*(x^2)
R4_growth <- VAR4/(VAR4 + Vresid_growth)

"Growth - September"; mean(R4_growth); quantile(R4_growth, probs = c(0.025, 0.975))
## [1] "Growth - September"
## [1] 0.9912901
##      2.5%     97.5% 
## 0.9876578 0.9941299